%% Parameter Sweep of ODEs
% This is a parameter sweep study of a 2nd order ODE system.
% m*X'' + b*X' + k*X = 0
% We solve the ODE for a time span of 0 to 25 seconds, with initial conditions X(0) = 0 and X'(0) = 1. We sweep the parameters "b" and "k" and record the peak values of X for each condition. At the end, we plot a surface of the results.
% We speed up the time to run our simulation by the use of the parfor construct.

% Further developed by HA 28.11.2016
type odesystem.m

%
%

clear
close all

%
% Underdamped -> vibration b^2 < 4 m k

m=1;
b=0:.1:5;
k=(b.^2)/4;

plot(b,k);grid on;xlabel('damping b');ylabel('stiffness k');
title('k=(b^2)/4, up left: vibration')
shg

%
%% Initialize parameters

m     =         1;  % mass
bVals =  .1:.1 :5;  % damping values
kVals = 0.5:.05:5;  % stiffness values

%
%

hold on
[bGrid, kGrid] = meshgrid(bVals, kVals);
BK=[bGrid(:) kGrid(:)];
plot(BK(:,1),BK(:,2),'.')

%
%

bb=BK(:,1);kk=BK(:,2);
vibr=(bb.^2 < 4*kk);
novibr=~vibr;
byes=bb(vibr);
bno=bb(novibr);
kyes=kk(vibr);
kno=kk(novibr);
hold on
plot(byes,kyes,'ok')
plot(bno,kno,'y.')

%
%

peakValsvibr=nan(length(byes),1);
peakValsno=nan(length(bno),1);

%
%% Let's make one call first for each case
% Note the handy handling of extra parameters.

[Tyes{1},Yyes{1}] = ode45(@(t,y) odesystem(t, y, m, byes(1), kyes(1)),...
    [0, 25], ...  % simulate for 25 seconds
    [0, 1]) ;      % initial conditions
[Tno{1},Yno{1}] = ode45(@(t,y) odesystem(t, y, m, bno(1), kno(1)),...
    [0, 25], ...  % simulate for 25 seconds
    [0, 1]) ;      % initial conditions

max(Yyes{1}(:,1))   % Max displacement, vibration case
min(Yyes{1}(:,2))   % Min velocity, vibration case
max(Yno{1}(:,1))    % Max displacement, no vibration case
min(Yno{1}(:,2))    % Min velocity, no vibration case

%
%% Parameter Sweep in serial

disp('Computing Parameter Sweep...');
% Preallocate
Tyes=cell(1,numel(byes));Yyes=cell(1,numel(byes));
tic;
for idx = 1:numel(byes)
  % Solve ODE
  [Tyes{idx},Yyes{idx}] = ode45(@(t,y) odesystem(t, y, m, byes(idx), ...
   kyes(idx)),[0, 25], ...  % simulate for 25 seconds
   [0, 1]) ;      % initial conditions
 peakValsvibr(idx)=max(Yyes{idx}(:,1));
end
toc   % 12.369371 Preallocation of cell-arrays above -> no effect

%
%

maxindyes=find(peakValsvibr==max(peakValsvibr))

%
%

tic
for idx = 1:numel(bno)
[Tno{idx},Yno{idx}] = ode45(@(t,y) odesystem(t, y, m, bno(idx), ...
   kno(idx)),[0, 25], ...  % simulate for 25 seconds
   [0, 1]) ;      % initial conditions

peakValsno(idx)=max(Yno{idx}(:,1));
end
toc   % 6.701297

%
%

maxindno=find(peakValsno==max(peakValsno))

%
%

subplot(2,1,1)
plot(Tyes{1},Yyes{1})
legend('y','y''')
title(['Vibration, max amplitude, b= ' num2str(byes(1))...
    ' k= ' num2str(kyes(1)) ])
subplot(2,1,2)
plot(Tno{1},Yno{1});
legend('y','y''')
title(['No vibration, max displacement, b= ' num2str(bno(1))...
    ' k= ' num2str(kno(1)) ])
shg
% Elapsed time is 13.047220 seconds

%
%

figure
hold on
for idx=[1 10 50 500]
    plot(Tyes{idx},Yyes{idx}(:,1))
end

%
%

figure
hold on
for idx=[1 10 50 500]
    plot(Tno{idx},Yno{idx}(:,1))
end

%
% Copyright 2013 The MathWorks, Inc.