
xini = [ 5.E+26 0];
h = 1/3;
ti = 0;
tf = 8;

% solve using Euler (eulsys.m)
[t, xeul] = eulsys(@decay,[ti tf],xini,h);

% solve using ODE45
[t, xode45] = ode45(@decay,[ti:h:tf],xini);

% solution for St-92
str92 = @(t,x) xini(1)*exp(-0.256*t); 

% plot
figure
plot(t,xeul,t,xode45')
hold on
plot(t,str92(t),'b.')
legend('St-92 (Eul)','Ytt-92 (Eul)','Str-92 (ODE45)','Ytt-92 (ODE45)','St-92 (analytical)')
xlabel('time (hours)')

% table
format short e
't St-92(Eul) Str-92(ODE45) St-92(analytical) rel.error.Eul(%) rel.err.ODE45(%)'
[t xeul(1,:)' xode45(:,1) str92(t) abs(xeul(1,:)'-str92(t))./str92(t) abs(xode45(:,1)-str92(t))./str92(t)]


