function [t,x,epsilon] = RK45adaptiveExercise5at10_3(int,x1,t1,t0,n)
%Input: Interval [a,b], "final" condition x1 at time t1, 
%   with t0 the time when the solution is to be computed (backwards),
%   and the number of steps n
%Output: the approximate solution x, at time T
%Example usage: [t,x,epsilon]=RK45adaptiveExercise5at10_3([1/2 3],6,3,1/2,100);
itmax=1.0e3;
epsmin=1.0e-10; 
epsmax=1.0e-8;
hmin=1.0e-4;
hmax=6.5e-2;
% To use the RKAdaptive45, in order to march backwards, one needs to switch
% the ends of the time interval, from [t0 t1] to [t1 t0]
[t,x,epsilon] = RKAdaptive45('(3*x)/t+4.5*t-13',[t1 t0],x1,n,itmax,epsmin,epsmax,hmin,hmax);

dt=1.e-3;
t2=2*t1;
yexact=exactsolution(t0:dt:t1);
yexact1=exactsolution(t0:dt:2*t1);
%

[t45b,y45b]= ode45(@rhs,[t1 t0],exactsolution(t1));%backward
[t15sb,y15sb]= ode15s(@rhs,[t1 t0],exactsolution(t1));%backward
[t23sb,y23sb]= ode23s(@rhs,[t1 t0],exactsolution(t1));%backward
[t45,y45]= ode45(@rhs,[t0 t1],exactsolution(t0));%forward
[t15s,y15s]= ode15s(@rhs,[t0 t1],exactsolution(t0));%forward
[t23s,y23s]= ode23s(@rhs,[t0 t1],exactsolution(t0));%forward
[ count, iter, xp, yp ] = ode_midpoint_VS ('3*(y/x)+(9/2)*x-13', t0, t1, exactsolution(t0), 1.e-4 , 1.e-6 );
%[t,x,epsilon] = RKAdaptive45('3*x/t+9/2*t.^2-13',[t0 t1],2.25,n,itmax,epsmin,epsmax,hmin,hmax);

figure(2)
hold all
plot(t45b,y45b,'r*',...
 t15sb,y15sb,'b*',...
 t23sb,y23sb,'g*',...
t0:dt:t1,yexact,'r-')
legend('ode45 backward','ode15s backward','ode23s backward','exact solution','Location','NorthEastOutside')
%
figure(3)
hold all
plot(t45,y45,'ro',...
t15s,y15s,'bo',...
t23s,y23s,'go',...
xp,yp,'y*-',...
t0:dt:t1,yexact,'r-')
legend('ode45 forward','ode15s forward','ode23s forward','midpoint','exact solution','Location','NorthEastOutside')
%
function zz = rhs(t,x)
zz = 3*(x/t)+(9/2)*t-13;
%
function z = exactsolution(t)
z = t.^3-(9/2)*t.^2+(13/2)*t;
