%
%   Program to solve the equation 
%
%   d^2y/dt^2 = -4y + sin(At)
%
%   This system is equivalent to 
%
%   dy/dt = z
%   dz/dt = -4y + sin(At)
%
global A;
                           

A  = input(' enter value of constant A: ');
y0 = input(' enter value of y(0): ');
z0 = input(' enter value of dy/dt(0), i.e. z(0): ');

t0 = 0;
tfinal = input(' enter final time: ');
delt = tfinal/100;
tspan=t0:delt:tfinal;

uin(1) = y0;                           % initial conditions
uin(2) = z0;                           %
[t,u] = ode45('order2_rhs',tspan,uin); % variable step RK4/5


y = u(:,1);                   % ode45 returns y in col.1 of u
z = u(:,2);                   % ode45 returns z in col.2 of u
%plot(t,y,t,yexact,'r+')
plot(t,y)
xlabel('time');
ylabel('y');
title('second order ODE ');