Task #1
Write a single function script to perform the Fourth Order Adams-Bashforth-Moulton algorithm on a system of ordinary differential equations (IVP). The program must be able to handle any number of differential equations (an IVP with any number of state variables).
Your function m-file will receive through its header the following problem parameters (in this
exact order):
a. a cell array containing the function handles of the equations describing the time derivatives of the several state variables
b. the initial augmented state-array – a numerical row array containing the needed starting values of the independent variable and the several state variables (in the same order as the function handles in the cell array from the first input argument)
c. the step size to be used, and
d. the number of steps to take

Function Output
Your file should not return anything to the calling program. Instead, it should
i. print to screen (the command window) the final state of the system, including t, and
ii. plot a trace for each of the dependent variables as a function of time.
Plot each trace on a separate graph. To facilitate plotting, store the results of all time steps in the augmented state-array
a. each row represents the system at a different point in time
b. column 1 contains the time, and each other column contains the temporal history of a particular state variable
Store only the final result of each time step, not the predicted or non-converged corrected values. Title each plot ‘Var1’, ‘Var2’, etc. Try help plot and help title in the MATLAB command window to learn more about plotting. num2str will also be useful for creating the labels themselves.

Hints for your Function Script
Handling a System of Un-predetermined Size
Technical
Because you will not know in advance the order of the IVP your function will be asked to solve, some of the system parameters will need to be passed in as arrays, rather than as scalars. The step size and number of steps will still be scalars, but there will be an un-predetermined number of state variables, each of which will have an associated starting value and derivative expression. The starting value of each state variable is a scalar, so the collection of them can be stored in a normal array. On the other hand, your function script must be passed a separate function handle for the equation for the time derivative of each state variable. To store function handles, we must use a cell array. So then, you must be able to evaluate a function handle stored in a cell array. We did this in Lab #4, but I cover it again later in this document.
NOTE: Because your augmented state array will include time as the first value, the function stored in the nth index of the cell array is actually the derivative of the state variable in the (n+1)th index of the augmented state array.

Function handles with unknown number of inputs
If you remember doing this for Modified Secant on a system of equations, then you may skip this section.
Up to now, you’ve probably only used function handles for equations that had a predetermined number of inputs. For this program, each derivative equation may be a function of the independent variable and any or all of the state variables. This is only a problem because you don’t know now (as you’re writing the function script) how many or which of the state variables will appear in each of the derivative equations. Fortunately, MATLAB allows us to put arrays into the equation function evaluation.
For example, if I define a function handle:
d1_by_dt = @(t,a) a(1) * a(3) / (2 + t);
when you ask MATLAB to evaluate the mathematical expression, MATLAB will expect to receive a scalar as the 1st argument, and an array of at least three terms as the 2nd argument. It will use the first and third terms in the array, along with the scalar, to evaluate the expression.
The augmented state array passed into the function already has t as the first value, so let’s just stick with that. Then our example equation becomes:
d1_by_dt = @(a) a(2) * a(4) / (2 + a(1));
You, of course, will have no prior knowledge of which state variables are used in each derivative expression. Therefore, for each derivative evaluation, you will need to submit an array containing the correct values of the dependent variable and all of the state variables, in the same order as the initial values in the passed-in augmented state array. When grading your code, it will be my responsibility to ensure that all of the derivative equations use the correct variables.

Function handles in cell arrays
The function handles of the derivative expressions of the several state variables will be passed to your function script in a cell array. The syntax for evaluating a function handle that is stored in a cell array is as follows:
value = array_name{n}(input1, input2, … , inputQ)
• array_name is the name of the cell array in which the handles are stored
• n is the index into the cell array, in which is stored the expression for (d/dt)y n
• input1 thru inputQ are the inputs to the mathematical expression represented by the function handle – in our case that will simply be the augmented state array
Please note that those are braces around the index ‘n’ and parentheses around the arguments.

The ABM4 process
To calculate the predicted value of each state variable, you are going to perform very similar math:
n+1 ^0y k = n y k + h((55/24)g^k n - (59/24)g^k n-1 + (37/24)g^k n-2 - (5/24)g^k n-3)

where the k superscript on each g term matches the state variable number (the k subscript on each y term). Take advantage of how MATALB deals with vectors to simplify the process of calculating the predictor for all the state variable and each subsequent corrector for all the state variables. Create an array that holds the g m evaluations for all of the state variables. You will have 3 or 4 of these, one for each timestep. Then you can …
temp_state = current_state + (h / 2) * slopes_array

Test your program
The spring-mass system from the Terminology section has been reproduced here.
Neglecting friction, the equations of motion of this system are:
ẍ1 = (-(k1+k2)x1 + k2x2)/m1
ẍ2 = (-(k2+k3)x2 + k2x1)/m2

where x1 and x2 are measured from the static equilibrium locations of the respective masses. To handle the 2nd order nature of the time derivatives, we must include additional state variables v1 and v2. If we arrange the state variables in the order x1, v1, x2, v2, the augmented state array becomes t, x1, v1, x2, v2, and the corresponding differential equations becomes:
dx1/dt = g1(a) = a(3)
dv1/dt = g2(a) = (-(k1+k2)a(2) + k2a(4))/m1
dx2/dt = g3(a) = a(5)
dv2/dt = g4(a) = (-(k2+k3)a(4) + k2a(2))/m2

Remember, a(1) is time, which does not appear in our equations of motion.
Create a main script (program) that performs the following tasks:
iii. define the system parameters (masses and spring stiffnesses)
iv. define the IVP functions above and store them in a cell array
v. set the initial condition of the augmented state array (time and all four state variables)
vi. call your RK4 function script
Set the system parameters to m1 = m2 = 1(kg) and k1 = k2 = k3 = 13.16(N/m). Set the initial
conditions to t = 0, x1(0) = 0.01(m), v1(0) = 0, x2(0) = -0.01(m), and v2(0) = 0. Call your RK4 function script with a time step of 0.001 seconds and 2000 time steps. Inspect the plots of x1 and x2. They should be mirrored sinusoidal functions with amplitudes of 0.01 and frequencies of 1Hz (2π rad/sec).

Analyze Your Custom Vibration System
You are going to customize the parameters of the spring mass system based upon your student ID number. Consider digits 4-5 as a two digit number X. Let digits 6-7 be a two digit number Y. If either X or Y are zero, pretend that they are 35.

Task #2: Apply Your Program
You will use your function to analyze a customized version of the vibration system depicted on the previous page. The system will be customized based on the last four digits of your student ID. Your goal will be to determine the second natural frequency and corresponding mode shape. For more details, refer to the Terminology section of this document.
Set the system parameters to m1 = 1(kg), m2 = (0.75-(X/200))(kg), k1 = k2 = 13.16(N/m), and k3 = (13.66+(Y/200))(N/m). Set the initial conditions to t = 0, x1(0) = -0.01(m), v1(0) = 0, and v2(0) = 0. Please note that the initial value of x1 is negative. Now, find the second natural frequency and second mode shape of your system by running your RK4 function multiple times, with different initial values for x2. I think the correct value should be > 0.01.
When you find the correct value of x2(0) each mass should demonstrate smooth sinusoidal motion. Inspect your table of state variables to determine the period of the oscillation. The period is the time between consecutive maximum displacements. The frequency in Hz is the inverse of the period. The mode shape is the vector [x1(0); x2(0)] scaled to have a Euclidean norm of unity (1).
Include in your Lab Report the following:
• the last four digits of your student ID number
• the values you used for X and Y
• the second natural frequency, in Hz, of your customized spring-mass system
• the normalized second mode shape of your customized spring-mass system
• plots of x1 and x2 as functions of time for two seconds
Provide the natural frequency and mode shape values to 5 significant figures.

Here is the RK4 chunk that is supposed to go on the top of the code for the function:
% In the code below:
%
% g is cell array containing the various dy/dt functions
%
% s_a is the augmented state array. It starts out as a row vector. It is
% passed in with the starting values of t and each state variable, in the
% same order as the function handles in the g cell array. It will
% eventually contain all of the values of t and the state variables over
% time.
%
% h is the step size

N = length(g);
slopes = zeros(4,N)

for iter = 1:2
state = s_a(iter,:);
for c = 1:N
slopes(1,c) = g{c}(state);
end

state = s_a(iter,:) + (h/2)[1, slopes(1,:)];
for c = 1:N
slopes(2,c) = g{c}(state);
end

state = s_a(iter,:) + (h/2)[1, slopes(2,:)];
for c = 1:N
slopes(3,c) = g{c}(state);
end

state = s_a(iter,:) + h*[1, slopes(3,:)];
for c = 1:N
slopes(4,c) = g{c}(state);
end

RK = (slopes(1,:) + 2slopes(2,:) + 2slopes(3,:) + slopes(4,:)) / 6;

s_a(iter+1,:) = s_a(iter,:) + h*[1, RK]

end

When this chunk is done, we should have three rows in the augmented state array, corresponding to our first three time points. It shouldn’t be changed in any way while making the function script.