Fix this MatLab code. There's an unrecognized function or variable "max_theta".


function IntegratingMethodsSHO(max_theta,method_name,dt_size)
% Solve Eq. of motion using a variety of integration methods
% ------------------------------------------------------------------------
% INPUT
% Maximum theta (Amplitude)                                   % in degrees 
% Method type: Euler
%              Euler-Cromer 
%              Runge-Kutta 
%              Verlet
%              leapfrog
% dt_size                                                     % in seconds
% 
% PROCESS
% Solve simple pendulum problem in the small angle approximation
% Use different type of numerical methods as specified by the user.
%
% methods = Euler, Euler-Cromer, Runge-Kutta, Verlet and leapfrog. 
%
% OUTPUT
% for selected method:
% Plot theta versus time for the selected method on figure 1. 
% Plot energy versus time for the selected method on figure 2. 
%
% for selected method = 'stable'
% Plot theta versus time for all methods except Euler on same figure 1. 
% Plot energy vs. time for all methods except for Euler on same figure 2.
%
% ------------------------------------------------ define fixed parameters 
  L = 2;                                % length of pendulum in units in m
  g = 10;                           % acceleration due to gravity in m/s^2
  max_theta = max_theta*pi/180;          % convert from degrees to radians 
% ----------------------------------------------------- input error checks
theta_max = 30*pi/180;                % maximum theta amplitude in radians
   if( nargin == 0 )
   max_theta = theta_max;                           % default = 30 degrees
   method_name = 'Euler-Cromer';                          % default method
   dt_size = 0.1;                                 % default time increment 
   elseif( nargin == 1 ) 
   method_name = 'Euler-Cromer';                          % default method
   dt_size = 0.1;                                 % default time increment
   elseif( nargin == 2 ) 
   dt_size = 0.1;                                 % default time increment
   end
% -------------------------------------------------------- check amplitude 
   if( isnumeric(max_theta) == 0 )
   disp('USAGE: IntegratingMethodsSHO(max_theta,method_name,dt_size)');
   disp('       max_theta must be numeric');
   error('ERROR: max_theta is not a numeric variable.');
   end
   if( max_theta < 0.0 )
   error('ERROR: Amplitude of oscillation cannot be less than 0');
   end
   if( max_theta > theta_max )
   error('ERROR: Amplitude of oscillation exceeds small angle approx');
   end
% --------------------------------------------------- check time increment
   if( isnumeric(dt_size) == 0 )
   disp('USAGE: IntegratingMethodsSHO(max_theta,method_name,dt_size)');
   disp('       dt_size must be numeric');
   error('ERROR: dt_size is not a numeric variable.');
   end
   if( dt_size < 1.0e-5 )
   error('ERROR: dt_size should be equal to or greater than 0.00001 sec');
   end
   if( dt_size > 1 )
   error('ERROR: dt_size should be equal to or less than 1 sec');
   end
% ---------------------------------------------------------- specification
spec = struct;
spec.L = L;
spec.g = g;
spec.m = 5;                                                  % units of kg
spec.theta0 = max_theta;
spec.omega0 = 0;
spec.dt = dt_size;
% ------------------------------------------------------ check method type
  switch lower(method_name)             % convert string to all lower case
    case 'euler'   % ---------------------------------------- Euler method
    [theta,energy,t] = GetInfoEuler(spec);
    figure; plot(t,theta );
    title('theta versus time graph');
    ylabel('theta (radians)');
    xlabel('time (seconds)');
    figure; plot(t,energy);
    title('Energy versus time graph');
    ylabel('energy (Joules)');
    xlabel('time (seconds)');
    case 'euler-cromer'  % --------------------------- Euler-Cromer method
    [theta,energy,t] = GetInfoEulerCromer(spec);
    figure; plot(t,theta );
    title('theta versus time graph');
    ylabel('theta (radians)');
    xlabel('time (seconds)');
    figure; plot(t,energy);
    title('Energy versus time graph');
    ylabel('energy (Joules)');
    xlabel('time (seconds)');
    case 'runge-kutta'  % ----------------------------- Runge-Kutta method
    [theta,energy,t] = GetInfoRungeKutta(spec);
    figure; plot(t,theta );
    title('theta versus time graph');
    ylabel('theta (radians)');
    xlabel('time (seconds)');
    figure; plot(t,energy);
    title('Energy versus time graph');
    ylabel('energy (Joules)');
    xlabel('time (seconds)');
    case 'verlet'  % --------------------------------------- verlet method
    [theta,energy,t] = GetInfoVerlet(spec);
    figure; plot(t,theta );
    title('theta versus time graph');
    ylabel('theta (radians)');
    xlabel('time (seconds)');
    figure; plot(t,energy);
    title('Energy versus time graph');
    ylabel('energy (Joules)');
    xlabel('time (seconds)');
    case 'leapfrog'  % ----------------------------------- leapfrog method
    [theta,energy,t] = GetInfoLeapfrog(spec);
    figure; plot(t,theta );
    title('theta versus time graph');
    ylabel('theta (radians)');
    xlabel('time (seconds)');
    figure; plot(t,energy);
    title('Energy versus time graph');
    ylabel('energy (Joules)');
    xlabel('time (seconds)');
    case 'stable'   % ---------------- plot results for all stable methods
    [theta1,energy1,t1] = GetInfoEulerCromer(spec);
    [theta2,energy2,t2] = GetInfoRungeKutta(spec);
    [theta3,energy3,t3] = GetInfoVerlet(spec);
    [theta4,energy4,t4] = GetInfoLeapfrog(spec);
    % ------------------------------------------------ plot theta vs. time
    figure;
    hold on;
    plot(t1,theta1,'b');
    plot(t2,theta2,'r');
    plot(t3,theta3,'k');
    plot(t4,theta4,'m');
    title('theta versus time graph');
    ylabel('theta (radians)');
    xlabel('time (seconds)');
    legend('Euler-Cromer','Runge-Kutta','Verlet','Leapfrog', ...
           'Location','southwest')
    hold off;
    % ----------------------------------------------- plot energy vs. time
    figure;
    hold on;
    plot(t1,energy1,'b');
    plot(t2,energy2,'r');
    plot(t3,energy3,'k');
    plot(t4,energy4,'m');
    title('Energy versus time graph');
    ylabel('energy (Joules)');
    xlabel('time (seconds)');
    legend('Euler-Cromer','Runge-Kutta','Verlet','Leapfrog', ...
           'Location','northwest')
    hold off;
    otherwise
    disp('OPTIONS: Euler, Euler-Cromer, Runge-Kutta, Verlet, Leapfrog');
    error('ERROR: Please select an implemented method');
  end 
end

function [theta,energy,t] = GetInfoLeapfrog(spec)
% Leapfrog method to calculate SHO of simple pendulum
% ------------------------------------------------- set time scale & steps
dt = spec.dt;
T = 2*pi*sqrt(spec.L/spec.g);
nsteps = floor(10*T/dt + 0.5);
if( mod(nsteps, 2) == 1 ) nsteps = nsteps + 1; end
% ----------------------------------------------------- preallocate memory
t = zeros(nsteps,1);
theta = t;
omega = theta;
energy = theta;
% ------------------------------------------------- set initial conditions
krestore = spec.g/spec.L; 
const = 0.5*spec.m*(spec.L*spec.g);
theta(1) = spec.theta0;                                 % at odd time step
omega(1) = spec.omega0;                               % In sync with theta
Econst = 0.5*spec.m*spec.L^2;
energy(1) = Econst*( omega(1)^2 + krestore*theta(1)^2 );
%   Apply Runge-Kutta method to calculate the second step for pendulum SHO
% ------------------------------------- get middle points (see A.16 in CP)
theta_m = theta(1) + 0.5*omega(1)*dt;      
omega_m = omega(1) - 0.5*krestore*theta(1)*dt;
% ---------------------------------------- add increments (see A.15 in CP)  
omega(2) = omega(1) - krestore*theta_m*dt;                % even time steps
theta(2) = theta(1) + omega_m*dt; 
t(2) = t(1) + dt;
energy(2) = Econst*( omega(2)^2 + krestore*theta(2)^2 );
dt2 = 2*dt;                                          % twice the time step
% ------------------------------------------ integrate equations of motion
nf = nsteps - 1;                               % one step less than needed
   for k=3:2:nf                                    % skip every other line
% -------------------------------------------------------- apply leap frog 
   theta(k) = theta(k-2) + omega(k-1)*dt2;                % k => odd steps   
   omega(k+1) = omega(k-1) - krestore*theta(k)*dt2;  
   t(k) = t(k-2) + dt2;
   omega_m = 0.5*( omega(k+1) + omega(k-1) );            % sync with theta
   energy(k) = Econst*( omega_m^2 + krestore*theta(k)^2 );
   end
% ----------------------------------------------- get energy on even steps
n2 = nsteps - 2;
   for k=4:2:n2
   theta(k) = 0.5*( theta(k+1) + theta(k-1) );
   energy(k) = Econst*( omega(k)^2 + krestore*theta(k)^2 );
   t(k) = t(k-2) + dt2;
   end
theta_m = theta(nf) + omega(nsteps)*dt; 
energy(nsteps) = Econst*( omega(nsteps)^2 + krestore*theta_m^2 );
t(nsteps) = t(n2) + dt2;
end

function [theta,energy,t] = GetInfoVerlet(spec)
% Verlet method to calculate SHO of simple pendulum
% ------------------------------------------------- set time scale & steps
dt = spec.dt;
T = 2*pi*sqrt(spec.L/spec.g);
nsteps = floor(10*T/dt + 0.5);
% ----------------------------------------------------- preallocate memory
t = zeros(nsteps,1);
theta = t;
omega = theta;
energy = theta;
% ------------------------------------------------- set initial conditions
krestore = spec.g/spec.L; 
const = 0.5*spec.m*(spec.L*spec.g);
theta(1) = spec.theta0;
omega(1) = spec.omega0;
energy(1) = const*theta(1)*theta(1);
Econst = 0.5*spec.m*spec.L^2;
energy(1) = Econst*( omega(1)^2 + krestore*theta(1)^2 );
%   Apply Runge-Kutta method to calculate the second step for pendulum SHO
% ------------------------------------- get middle points (see A.16 in CP)
theta_m = theta(1) + 0.5*omega(1)*dt;      
omega_m = omega(1) - 0.5*krestore*theta(1)*dt;
% ---------------------------------------- add increments (see A.15 in CP)  
omega(2) = omega(1) - krestore*theta_m*dt;  
theta(2) = theta(1) + omega_m*dt;   
t(2) = t(1) + dt; 
energy(2) = Econst*( omega(2)^2 + krestore*theta(2)^2 );
Adt2 = krestore*dt*dt; 
divide_by_2dt = 1.0/(2*dt);
% ------------------------------------------ integrate equations of motion 
   for k=3:nsteps
% ---------------------------------- apply 2nd order two step relationship
   theta(k) = ( 2*theta(k-1) - theta(k-2) ) - Adt2*theta(k-1);   
   omega(k-1) = ( theta(k) - theta(k-2) )*divide_by_2dt;  
   t(k) = t(k-1) + dt; 
   energy(k-1) = Econst*( omega(k-1)^2 + krestore*theta(k-1)^2 );  
   end
% --------------------------- apply Runge-Kutta estimate for omega(nsteps)
theta_m = theta(nsteps-1) + 0.5*omega(nsteps-1)*dt;
omega(nsteps) = omega(nsteps-1) - krestore*theta_m*dt;
energy(nsteps) = Econst*( omega(nsteps)^2 + krestore*theta(nsteps)^2 );
end

function [theta,energy,t] = GetInfoRungeKutta(spec)
% Runge-Kutta method to calculate SHO of simple pendulum
% ------------------------------------------------- set time scale & steps
dt = spec.dt;
T = 2*pi*sqrt(spec.L/spec.g);
nsteps = floor(10*T/dt + 0.5);
% ----------------------------------------------------- preallocate memory
t = zeros(nsteps,1);
theta = t;
omega = theta;
energy = theta;
% ------------------------------------------------- set initial conditions
krestore = spec.g/spec.L; 
const = 0.5*spec.m*(spec.L*spec.g);
theta(1) = spec.theta0;
omega(1) = spec.omega0;
Econst = 0.5*spec.m*spec.L^2;
energy(1) = Econst*( omega(1)^2 + krestore*theta(1)^2 );
% ------------------------------------------ integrate equations of motion
   for k=2:nsteps
% ------------------------------------- get middle points (see A.16 in CP)
   theta_m = theta(k-1) + 0.5*omega(k-1)*dt;      
   omega_m = omega(k-1) - 0.5*krestore*theta(k-1)*dt;
% ---------------------------------------- add increments (see A.15 in CP)  
   omega(k) = omega(k-1) - krestore*theta_m*dt;  
   theta(k) = theta(k-1) + omega_m*dt;   
   t(k) = t(k-1) + dt; 
   energy(k) = Econst*( omega(k)^2 + krestore*theta(k)^2 );  
   end
end

function [theta,energy,t] = GetInfoEulerCromer(spec)
% Euler-Cromer method to calculate SHO of simple pendulum
% ------------------------------------------------- set time scale & steps
dt = spec.dt;
T = 2*pi*sqrt(spec.L/spec.g);
nsteps = floor(10*T/dt + 0.5);
% ----------------------------------------------------- preallocate memory
t = zeros(nsteps,1);
theta = t;
omega = theta;
energy = theta;
% ------------------------------------------------- set initial conditions
krestore = spec.g/spec.L; 
const = 0.5*spec.m*(spec.L*spec.g);
theta(1) = spec.theta0;
omega(1) = spec.omega0;
Econst = 0.5*spec.m*spec.L^2;
energy(1) = Econst*( omega(1)^2 + krestore*theta(1)^2 );
% ------------------------------------------ integrate equations of motion
   for k=2:nsteps
   omega(k) = omega(k-1) - krestore*theta(k-1)*dt;
   theta(k) = theta(k-1) + omega(k)*dt;            % see example 3.2 in CP
   t(k) = t(k-1) + dt; 
   energy(k) = Econst*( omega(k)^2 + krestore*theta(k)^2 );
   end
end

function [theta,energy,t] = GetInfoEuler(spec)
% function trajectory = GetInfoEuler(spec)
% Euler method to calculate SHO of simple pendulum
% ------------------------------------------------- set time scale & steps
dt = spec.dt;
T = 2*pi*sqrt(spec.L/spec.g);
nsteps = floor(10*T/dt + 0.5);
% ----------------------------------------------------- preallocate memory
t = zeros(nsteps,1);
theta = t;
omega = theta;
energy = theta;
% ------------------------------------------------- set initial conditions
krestore = spec.g/spec.L; 
const = 0.5*spec.m*(spec.L*spec.g);
theta(1) = spec.theta0;
omega(1) = spec.omega0;
Econst = 0.5*spec.m*spec.L^2;
energy(1) = Econst*( omega(1)^2 + krestore*theta(1)^2 );
% ------------------------------------------ integrate equations of motion
   for k=2:nsteps
   omega(k) = omega(k-1) - krestore*theta(k-1)*dt;
   theta(k) = theta(k-1) + omega(k-1)*dt;          % see example 3.1 in CP
   t(k) = t(k-1) + dt; 
   energy(k) = Econst*( omega(k)^2 + krestore*theta(k)^2 );  
   end
end