%% Detta program löser mu'' + cu' + ku = Asin(oemga*t)  omskrivet
%% på standardform u'' = f(t,u,u')= 1/m*(Asin(oemga*t) -  cu' - ku)
%% med begynnelsevärden u(0)=0 och u'(0)=0 mha Eulers metod samt ODE45


clear all

global m c k A omega    % Använd globala variabler för att kunna 
                        % använda dom inne i en funktion 
                        % se fun_ex.m

 
t0=0;                   % Begynnelsedata
w0=[0; 0];

m=10;                   % Giva parameterar
c=0.1;
k=0.8;
A=4;
omega=0.5;

tend = 200;            % Sluttid
h=0.1;                 % Steglängd

w=w0;
t=t0;
tplot=t';              % Vektorer för att spara lösn och tidpunkter
wplot=w';

nosteps = tend/h;
% Eulers metod

for i=1:nosteps
  
  F=[w(2); 1/m*(A*sin(omega*t)-c*w(2)-k*w(1))];
  
  %% Alternativt kan man använda funktionen som vi skapat för ODE45
  %% istället för att räkna ut funktionen enl ovan
  %% F=fun_ex(t,w);
  
  w = w + h*F;
  t=t+h;
  tplot=[tplot;t'];
  wplot=[wplot;w'];
end

subplot(211)
plot(tplot,wplot(:,1),'b'), hold on 
subplot(212)
plot(tplot,wplot(:,2),'b'), hold on 


opt=odeset('RelTol',1e-6);
[tout,wout]=ode45('fun_ex',[0 tend],[0;0]);
subplot(211)
plot(tout,wout(:,1),'r--')
subplot(212)
plot(tout,wout(:,2),'r--')

