
function z=f(t,y)
  z = exp(t)+3*y
endfunction

function dz=df(t,y) // dans la méthode de Newton la pente est donnee par df/dy
  dz = 3
endfunction

function x=newton(t,s,yn)
  x=yn
  for j=1:10
    x=x-(yn+(s-t)*(f(t,yn)+f(s,x))/2-x)/( (s-t)/2*df(s,x)-1)
  end
endfunction

function y=eulerimp(t,y0)
  n=prod(size(t))
  y=zeros(n,1)
  y(1)=y0
  for i=1:n-1
    y(i+1)=newton(t(i),t(i+1),y(i))
  end
endfunction

eulerimp([0:1/8:1],12)

tps=[0:1/8:1];
plot(tps,eulerimp(tps,12));
plot(tps,-0.5*exp(tps)+12.5*exp(3*tps),'r-');

