funcprot(0)
// euler_exp impératif
function  E = euler_imp(f,a,b,N,yo)
  h = (b-a)/N
  tn = a
  un = yo
  X = tn
  Y = un
  for k = 1:N
    un = un + h*f(tn,un)
    tn = tn + h
    X = [X,tn]
    Y = [Y,un]
  end
  E = [X;Y]
endfunction


function plot_euler_imp(f,a,b,N,yo)
  e = euler_imp(f,a,b,N,yo)
  plot(e(1,:) , e(2,:))
endfunction

// euler_exp récursif
function  E = euler_rec_les_y(f,x,y,xmax,h,liste)
  if x > xmax then E = liste
  else E = euler_rec_les_y(f,x+h,y+h*f(x,y),xmax,h,[liste,y+h*f(x,y)])
  end
endfunction

function euler_rec(f,a,b,N,yo)
  h = (b-a)/N
  x = [a:h:b]
  plot(x,euler_rec_les_y(f,a,yo,b,h,[yo]))
endfunction

//plot_euler_imp(deff('[u] = f(x,y)', 'u = y') ,0,2,20,1),plot([0:0.1:2],exp([0:0.1:2]),":r")

//euler_rec(deff('[u] = f(x,y)', 'u = y') ,0,2,20,1),plot([0:0.1:2],exp([0:0.1:2]),":r")


//max(abs(euler_rec_les_y(deff('[u] = f(x,y)', 'u = y'),0,1,2,0.01,[1]) - exp([0:0.01:2])))

//e = euler_imp(deff('[u] = f(x,y)', 'u = y'),0,2,4096,1);
//max(abs((e(2,:) - exp(e(1,:)))))



// euler2_exp impératif
function  E = euler2_imp(f,a,b,N,yo)
  h = (b-a)/N
  tn = a
  un = yo
  X = tn
  Y = un
  for k = 1:N do
    un = un + h*f(tn + 0.5*h , un + 0.5*h*f(tn,un) )
    tn = tn + h
    X = [X,tn]
    Y = [Y,un]
  end
  E = [X;Y]
endfunction

function plot_euler2_imp(f,a,b,N,yo)
  e = euler2_imp(f,a,b,N,yo)
  plot(e(1,:) , e(2,:))
endfunction

// euler2_exp récursif
function  E = euler2_rec_les_y(f,x,y,xmax,h,liste)
  if x > xmax then E = liste
  else E = euler2_rec_les_y(f,x+h,y+h*f(x+0.5*h,y+0.5*h*f(x,y)),xmax,h,[liste,y+h*f(x+0.5*h,y+0.5*h*f(x,y))])
  end
endfunction

function euler2_rec(f,a,b,N,yo)
  h = (b-a)/N
  x = [a:h:b]
  plot(x,euler2_rec_les_y(f,a,yo,b,h,[yo]))
endfunction


//plot_euler2_imp(deff('[u] = f(x,y)', 'u = y') ,0,2,20,1),plot([0:0.1:2],exp([0:0.1:2]),":r")

//euler2_rec(deff('[u] = f(x,y)', 'u = y') ,0,2,20,1),plot([0:0.1:2],exp([0:0.1:2]),":r")


//max(abs(euler2_rec_les_y(deff('[u] = f(x,y)', 'u = y'),0,1,2,0.01,[1]) - exp([0:0.01:2])))
//e1 = euler_imp(deff('[u] = f(x,y)', 'u = y'),0,2,4096,1);
//e2 = euler2_imp(deff('[u] = f(x,y)', 'u = y'),0,2,4096,1);
//max(abs((e1(2,:) - exp(e1(1,:)))))
//max(abs((e2(2,:) - exp(e2(1,:)))))


//RK4


function R = rk4(f,a,b,N,yo)
  h = (b-a)/N
  tn = a
  un = yo
  X = tn
  Y = un
  for i = 1:N
    k1 = h*f(tn , un)
    k2 = h*f(tn+0.5*h , un+k1/2)
    k3 = h*f(tn+0.5*h , un+k2/2)
    k4 = h*f(tn+h , un+k3)
    tn = tn + h
    un = un + (k1 + 2*k2 + 2*k3 + k4)/6
    X = [X,tn]
    Y = [Y,un]
  end
  R = [X;Y]
endfunction

deff('[f] = f(x,y)', 'f = -150*y + 30');

//e = euler_imp(f,0,1,50,1/5+1D-15);
//e2 = euler2_imp(f,0,1,50,1/5+1D-15);
//r = rk4(f,0,1,50,1/5+1D-15);
//plot(r(1,:),r(2,:))
//max(abs((r(2,:) - exp(r(1,:)))))


//t = 0:0.02:1;
//uo = 1/5+1D-15;
//to = 0;
//u = ode(uo,to,t,f);
//plot(t,u)


// systèmes et ordre 2
// y1' = 3y1 -4y2
// y'2 = y1 + 3y2
// y1(0) = 2  y2(0) = 0


//deff('dy=f(t,y)' , 'dy(1) = 3*y(1) - 4*y(2) ; dy(2) = y(1) + 3*y(2)');
//t0 = 0 ; y0 = [2 ; 0] ; t = 0: 0.1 : 2 ;
//z = ode(y0,t0,t,f) ;
//plot(t,z(1,:),t,z(2,:))
//plot(t,2*cos(2*t) .* exp(3*t),"ro")

//pendule
deff('dy=f(t,y)' , 'dy(1) = y(2) ; dy(2) = -y(2) -sin(y(1))');
t0 = 0 ; y0 = [0 ; 1] ; t = 0: 0.1 : 10 ;

z = ode(y0,t0,t,f) ;
plot(t,z(1,:),t,z(2,:))


// portrait de phase
plot(z(1,:),z(2,:))
