clear all
clf
close all
clc

alpha_vector=linspace(0,15,100);
S_vector=zeros(1,100);

for j=1:length(alpha_vector)
    
    alpha=alpha_vector(j);
    t_end=1;
    
    dt = 1e-3;
    y0 = [0; alpha];
    y_old=y0;
    y_sol=zeros(2,t_end/dt+1);
    y_sol(:,1)=y0;
    
    for i=1:t_end/dt
        
        k1=dt.*RHS(y_old);
        k2=dt.*RHS(y_old+k1);
        y_new = y_old + (k1+k2)./2;
        
        y_sol(:,i+1)=y_new;
        
        y_old=y_new;
        
    end
    
    x=linspace(0,t_end,t_end/dt+1);
    y_exact=alpha.*sin(x);
    
%     figure(2)
%     hold on
%     plot(x,y_sol(1,:),'k')
%     plot(x,y_exact,'--m')
%     plot(x,y_sol(2,:),'r')
%     grid on
%     box on
    
    S=y_sol(2,end)-6;
    error=norm(y_sol(1,:)-y_exact,2);
    disp(['The error is      ' num2str(error)])
    disp(['The value of S is ' num2str(S)])
    
    S_vector(j)=S;
    
end

figure(1)
hold on
plot(alpha_vector,S_vector,'.k')
grid on
box on
xlabel('\alpha')
ylabel('S')