%% Solver Selection - Fixed Step

% Here is the most efficient way to choose the best fixed-step solver for
% your model experimentally. 
%%
% 1.	First, use one of the variable-step solvers
% to simulate your model to the level of accuracy that you desire. This
% will give you an idea of what the simulation results should be.
opts = simset('Solver','ode45','MaxStep','auto','RelTol',1e-3,'AbsTol',1e-3);
sim('sldemo_absbrake_digest4p1',[],opts)

%%
% Compare the results to ODE23 with the same tollerances.

opts = simset('Solver','ode23','MaxStep','auto','RelTol',1e-3,'AbsTol',1e-3);
sim('sldemo_absbrake_digest4p1',[],opts)

% 2. Next, use ode1 to simulate your model at the default step size for your
% model. Compare the results of simulating your model with ode1 with the
% results of simulating with the variable-step solver. If the results are
% the same within the specified level of accuracy, you have found the best
% fixed-step solver for your model, namely ode1. That's because ode1 is the
% simplest of the Simulink fixed-step solvers and hence yields the shorted
% simulation time for the current step size.
%
% If ode1 does not give accurate results, repeat the preceding steps with
% the other fixed-step solvers until you find the one that gives accurate
% results with the least computational effort. The most efficient way to do
% this is to use a binary search technique. First, try ode3. If it gives
% accurate results, try ode2. If ode2 gives accurate results, it is the
% best solver for your model; otherwise, ode3 is the best. If ode3 does not
% give accurate results, try ode5. If ode5 gives accurate results, try
% ode4. If ode4 gives accurate results, select it as the solver for your
% model; otherwise, select ode5. 
%
% If ode5 does not give accurate results, reduce the simulation step size
% and repeat the preceding process. Continue in this way until you find a
% solver that solves your model accurately with the least computational
% effort.

%% load up the variable step solver data

load absbrake_variable_step_solver

%% Plot all solvers time trace side by side

figure('Name','ABS Solver Step Size','NumberTitle','off');
solvers = {'ode45';'ode23';'ode113';'ode15s';'ode23s';'ode23t';'ode23tb'};
for i = 1:7
    subplot(4,2,i)
    plot(data.Ww.(solvers{i}).Time(2:end),diff(data.Ww.(solvers{i}).Time))
    title([solvers{i} ' - ' num2str(length(data.Ww.(solvers{i}).Time)) ' Simulation Time Steps'])
    ylabel('Step Size(secs)')
    xlabel('Time(secs)')
end
%% Compare Ww for all solvers, overlay
figure('Name','ABS Wheel Speed','NumberTitle','off');
axes
solvers = {'ode45';'ode23';'ode113';'ode15s';'ode23s';'ode23t';'ode23tb'};
hold on
inputList = cell(1,length(solvers)*2);
for i = 1:length(solvers)
    n = 2*i-1;
    inputList{n} = data.Ww.(solvers{i}).Time;
    inputList{n+1} = data.Ww.(solvers{i}).Data;
end
plot(inputList{:})
legend(solvers{:})
title('Wheel Speeds From Different Solvers')
ylabel('Speed(rad/sec)');
xlabel('Time(secs)')
grid on

%% Compare Ww for a couple fixed step solvers, overlay
figure('Name','ABS Wheel Speed','NumberTitle','off');
axes
solvers = {'ode1';'ode2';'ode3';'ode4';'ode5'};
hold on
inputList = cell(1,length(solvers)*2);
for i = 1:length(solvers)
    n = 2*i-1;
    inputList{n} = data.Ww.(solvers{i}).Time;
    inputList{n+1} = data.Ww.(solvers{i}).Data;
end
plot(inputList{:})
legend(solvers{:})
title(['Wheel Speeds From Different Solvers' 10 'Fixed Step Ts = ' num2str(diff(inputList{1}(1:2)))])
ylabel('Speed(rad/sec)');
xlabel('Time(secs)')
grid on


%% Compare Wwdot for a couple fixed step solvers, overlay
figure('Name','ABS Wheel Accel','NumberTitle','off');
axes
solvers = {'ode1';'ode2';'ode3';'ode4';'ode5'};
hold on
inputList = cell(1,length(solvers)*2);
for i = 1:length(solvers)
    n = 2*i-1;
    inputList{n} = data.Wwdot.(solvers{i}).Time;
    inputList{n+1} = data.Wwdot.(solvers{i}).Data;
end
plot(inputList{:})
legend(solvers{:},'Location','SouthWest')
title(['Wheel Acceleration From Different Solvers' 10 'Fixed Step Ts = ' num2str(diff(inputList{1}(1:2)))])
ylabel('Acceleration(rad/sec^2)');
xlabel('Time(secs)')
grid on