function odeexp
% solve dy/dt = t^3 - y^3 using euler's method
dt = 0.0156;
t(1) = 0;
ye(1) = 0;
n = 9/dt;
for i = 2:n
     ye(i) = ye(i-1) + dt*(t(i-1)^3 - ye(i-1)^3);
     t(i) = t(i-1) + dt;
end;

% solve dy/dt = t^3 - y^3 using 2nd order RK method
t(1) = 0;
y(1) = 0;
yt(1) = 0;
n = 9/dt;
for i = 2:n
     k1   = t(i-1)^3 - yt(i-1)^3;
     yt(i)= yt(i-1) + dt*k1;
     t(i) = t(i-1) + dt;
     k2   = t(i)^3 - yt(i)^3;
     y(i) = y(i-1) + 0.5*(k1+k2)*dt;
end;

figure(1);clf;
h=plot(t,ye,t,y,'r--');
set(h,'Linewidth',2);
set(gca,'Fontsize',20);
legend('ye','yr',2);
xlabel('time');
ylabel('y');
axis([0 9 0 10]);
grid;

pause;

clear all;

% solve y' = y(1-y) using euler's method
dt = 0.0156;
t(1) = 0;
ye(1) = 0.1;
n = 9/dt;
for i = 2:n
     ye(i) = ye(i-1) + dt*ye(i-1)*(1-ye(i-1));
     t(i) = t(i-1) + dt;
end;

% solve y' = y(1-y) using 2nd order RK method
t(1) = 0;
y(1) = 0.1;
yt(1) = 0.1;
for i = 2:n
     k1   = yt(i-1)*(1-yt(i-1));
     yt(i)= yt(i-1) + dt*k1;
     t(i) = t(i-1) + dt;
     k2   = yt(i)*(1-yt(i));
     y(i) = y(i-1) + 0.5*(k1+k2)*dt;
end;

figure(2);clf;
h=plot(t,ye,t,y,'r--');
set(h,'Linewidth',2);
set(gca,'Fontsize',20);
legend('ye','yr',2);
xlabel('time');
ylabel('y');
grid;

pause;

% solve y' = y(1-y) using ode45
y0 = 0.1;
options = odeset('RelTol',1e-4,'AbsTol',1e-4);
[tsol,ysol] = ode45(@rhs1,[0 9],y0,options);

figure(2); hold on;
h=plot(tsol,ysol,'c');
set(h,'Linewidth',2);
set(gca,'Fontsize',20);
xlabel('time');
ylabel('y');

pause;

ta = [0:0.05:9];
fa = 1./(1+9*exp(-ta));

figure(2); hold on;
h=plot(ta,fa,'m');
set(h,'Linewidth',2);
set(gca,'Fontsize',20);
xlabel('time');
ylabel('y');

pause;

fantsol = interp1(ta,fa,tsol);
ytsol     = interp1(t,y,tsol);
yetsol    = interp1(t,ye,tsol);

figure(3);
h=plot(tsol,abs(fantsol-ysol),tsol,abs(fantsol-ytsol),...
     tsol,abs(fantsol-yetsol));
set(h,'Linewidth',2);
set(gca,'Fontsize',20);
legend('ode45','rk','euler');
xlabel('time');
ylabel('error');
grid;

pause;
clear all;

% solve y' = t^3 - y^3 using ode45
y0 = 0.0;
options = odeset('RelTol',1e-6,'AbsTol',1e-6);
[tsol,ysol] = ode45(@rhs2,[0 40],y0,options);

figure(1); hold on;
h=plot(tsol,ysol,'c');
set(h,'Linewidth',2);
set(gca,'Fontsize',20);
xlabel('time');
ylabel('y');


function F = rhs1(t,y)

F = y(1)*(1-y(1)); 


function F = rhs2(t,y)

F = t^3 - y(1)^3;