%Gerardo

function EulerMethod

%% Section 2.7: 

%% Example 1: Linear Equation

%dy/dx=3-y+t

t_0=0; y_0=-1;

t_end=1;

Nsol=10; dt = t_end/Nsol;

t_sol = zeros(Nsol+1,1);
y_sol = zeros(Nsol+1,1);
y_exact = zeros(Nsol+1,1);

t_sol(1,1) = t_0;
y_sol(1,1) = y_0;
y_exact(1,1) = y_0;
for it=1:Nsol
    y_n = y_sol(it,1);
    t_n = t_sol(it,1);
    
    t_np1 = t_n+dt;
    m = 3-y_n+t_n;
    y_np1 = y_n+m*(t_np1-t_n);
    y_sol(it+1,1) = y_np1;
    t_sol(it+1,1) = t_np1;
    
    y_exact(it+1,1) = 2+t_np1+(y_0-2)*exp(-t_np1);
end

y_sol(Nsol+1,1)
y_exact(Nsol+1,1)
y_exact(Nsol+1,1)-y_sol(Nsol+1,1)


clf
figure(1)
hold on
plot(t_sol,y_exact,'b -','LineWidth',2);
plot(t_sol,y_sol,'r *','LineWidth',1);
stop

%% Example 2: Exact equation

%dy/dx=3-y+t

N=100;

t_ini=0; t_end=10;
y_ini=-3; y_end=0;

dt=(t_end-t_ini)/N; 
dy=(y_end-y_ini)/N; 

f=zeros(N+1,N+1,2); t_axis=zeros(N+1,1); y_axis=zeros(N+1,1);
g=zeros(N+1,N+1);
for iy=1:N+1
    y = y_ini+dy*(iy-1);
    y_axis(iy,1) = y;
    for it=1:N+1
        t=t_ini+dt*(it-1);
        t_axis(it,1) = t;
        f(it,iy,1) = cos(t)+t*exp(y)-2*y;
        f(it,iy,2) = y*sin(t)-exp(y);
        
        g(it,iy) = y*cos(t)+t*exp(y)-y^2;
    end
end

Feq=5;
% clf
% figure(1)
hold on
quiver(t_axis(1:Feq:N+1),y_axis(1:Feq:N+1),transpose(f(1:Feq:N+1,1:Feq:N+1,1)),transpose(f(1:Feq:N+1,1:Feq:N+1,2)),1)
plot(t_axis,y_axis*0)
plot(t_axis*0,y_axis)
axis([t_ini t_end y_ini y_end])
axis on
contour(t_axis,y_axis,transpose(g),100)
%stop

%

t_0=0; y_0=-1;
Nsol=100; dt = t_end/Nsol;

t_sol = zeros(Nsol+1,1);
y_sol = zeros(Nsol+1,1);

t_sol(1,1) = t_0;
y_sol(1,1) = y_0;
for it=1:Nsol+1
    y_n = y_sol(it,1);
    t_n = t_sol(it,1);
    
    t_np1 = t_n+dt;
    m = (y_n*sin(t_n)-exp(y_n))/(cos(t_n)+t_n*exp(y_n)-2*y_n);
    y_np1 = y_n+m*(t_np1-t_n);
    y_sol(it+1,1) = y_np1;
    t_sol(it+1,1) = t_np1;
end

plot(t_sol,y_sol,'r -','LineWidth',2);

