#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define DIVN 1000
#define M 1
#define DX 30

void euler(double [], double [][M], double, int, int);
void euler_sub(double, double [], double []);

int main(void)
{
  double h, x[DIVN], y[DIVN][M];
  int i, n;

  printf("\n** y' - {e^(0.01/(x-0.5)^2+0.002)}+5e^y = 0 **\n");
  x[0] = y[0][0] = 0; n = DX; h = 1. / DX;
  printf("X の初期値 X0 = %f\n", x[0]);
  printf("Y の初期値 Y0 = %f\n", y[0][0]);
  printf("刻み値 N = %d\n", n);
  printf("進み幅 H = %f\n", h);

  euler(x, y, h, n, M);
  printf("\%15s%15s\n\n","X","Y");
  for (i = 0; i <= n; i++)
    printf("%17.8e%17.8e\n", x[i], y[i][0]);

  return 0;
}

/*
 * オイラー法
 */
void euler(double x[], double y[][M], double h, int n, int m)
{
  int i, k;
  double f[DIVN];
  
  for (k = 0; k <= n; k++) {
    x[k] = x[0] + k * h;
    euler_sub(x[k], y[k], f);
    for (i = 0; i < m; i++)
      y[k + 1][i] = y[k][i] + h * f[i];
  }
}

void euler_sub(double x, double y[], double f[])
{
  f[0] = exp(0.01 / pow(x - 0.5, 2) + 0.002) + 5 * exp(y[0]);
}

y'=e^{0.01/(x-0.5)^2 + 0.002} + 5e^{y(x)}                  -(1)
初期条件 y(0)=0                                            -(2)
近似式   (y(x+h)-y(x)) / h = ｛f(x,y) + f(x+h, y(x+h))｝/2 -(3)
刻み幅 ｈ＝1/30

(1)(3)において
y(x)           を y(n)    
y(x+h)         を y(n+1) 
f(x, y)        を f(n) 
f(x+h, y(x+h)) を f(n+1)   (n = 0,1,2 ...)
と置き換えると

y'(n) = f(n) = e^{0.01/(x(n)-0.5)^2 + 0.002}  +  5e^{y(n)} -(4)
(y(n+1) - y(n)) / h = (f(n) + f(n+1)) / 2  (h=x(n+1)-x(n)) -(5)
となる

(5)の近似式を(4)に適用すると
(y(n+1) - y(n)) / h =  (e^{0.01/(x(n)-0.5)^2 + 0.002} + 5e^{y(n)}
                        + e^{0.01/(x(n+1)-0.5)^2 + 0.002}  +  5e^{y(n+1)}) / 2  -(6)
となる

(6)を変型すると
0 = y(n+1)
    - y(n)
    - h/2( e^{0.01/(x(n)-0.5)^2   + 0.002}  
           +  5e^{y(n)} 
           + e^{0.01/(x(n+1)-0.5)^2 + 0.002} ) 
    - h/2 * 5e^{y(n+1)}
   = z(y(n+1))     -(7)
となる

x(n), y(n), h, x(n+1)=x(n)+h の値がわかっていれば
z(y(n+1)) = 0 を y(n+1) について解くことができる(ニュートン法を使う)

n = 0 のとき
x(n) = x(0) = 0
y(n) = y(0) = 0
h = 1/30
x(n+1) = x(1) = x(0) + h = h
これらの値を(7)に代入して x(1)=0+h のときの y(1) の値を求める

n = 1 のときも同様の操作を繰り返す


これで出来ないでしょうか?

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "chinougl.h"
#define EXP 2.718281828
#define SC 0.8

void mdeul(int n, double x0, double y0, double xn);
void gr(int i, double x, double y);
typedef struct{
    double fnc;
    double dfn;
}FV;

double c=1.0;
int main(){
int i,n=30;
double x0=0.0,y0=0.0,xn=1.0;
mdeul(n, x0, y0, xn);c++;
//for(i=0;i<n;n++){x[i]=x0;y[i]=y0;}
gotogl();
return 0;
}

FV fn(double x, double y, double z){
    FV ff;
    ff.fnc=z-pow(EXP,0.01/(pow(x-0.5,2)+0.002))+5.0*pow(EXP,y);
    ff.dfn=1;
    return ff;
}

double F(double x, double y){
    int n=0,nmax=30;
    double er,e=1.0e-6;
    static double z=0.0;
    static int k=0;
    FV func;
    do {
        func=fn(x,y,z);
        er=func.fnc/func.dfn;
        z=z-er;
        n++;
    }while(fabs(er)>e && n<nmax);
    if(n==nmax)
    gotogl();
k++;
return z;
}


void mdeul(int n, double x0, double y0, double xn){
        int i;
        double x,xp,y,yp,h;
x=x0;
y=y0;
h = (xn - x0) / n;
for(i=0;i<n;i++){
    yp = y + F(x, y) * h;
    xp = x + h;
    y += (h / 2) * (F(x, y) + F(xp, yp));
    x += h;
    gr(i,x*SC,y*SC);
  }
    return;
}

void gr(int i, double x, double y){

static double xo,yo;

if (i){
      line(xo/2-1*SC,yo/2*SC,x/2-1*SC,y/2*SC);
      switch((int)c){
   case 1: color(0,0,1);break;
   case 2: color(1,0,0);break;
   case 3: color(0,1,0);break;
   case 4: color(1,1,1);
   }
   }

xo=x;
yo=y;