I'm encountering a problem in my code. When I try to plot the functions illustrated in my code, I draw a blank screen from the function. The code is meant to solve a mass-spring-dapshot system using Forward, Backward, and Crank_Nicolson. My code is below:

% Using the Forward/Backward Euler method for % solving a mass-spring-dashpot system. clear; % clears all memory clc; % clears Command Window screen close all; % closes all open figure windows % 1.1 Initialize dt = input('Insert time step: '); % choose delta t t_start = 0; % starting time t_end = 50; % end time nstep = ( t_end - t_start ) / dt;% number of time steps x = zeros(2,1); % x(1): displacement, x(2): velocity b = zeros(2,1); % right-hand-side vector I = eye(2); % identity matrix all_t = ( 1 : nstep ) * dt; % all time points all_Euler_x = zeros( nstep , 2 ); % Euler solution - pre-allocation all_Backward_y = zeros( nstep, 2); %Backward Euler - pre-allocation all_CN_z = zeros( nstep, 2); % Crank-Nicolson - pre-allocation all_exact_x = zeros( nstep , 2 ); % exact solution - pre-allocation initial_displacement = 1.0; % initial displacement initial_velocity = 0.0; % initial velocity % 1.2 system characteristics m = 1.0; % mass c = 0.05; % damping coefficient k = 4.0; % stiffness A(1,1) = 0.0 ; A(1,2) = 1.0; A(2,1) = -k/m ; A(2,2) = -c/m; % 1.3 initial conditions x(1) = initial_displacement; % initial displacement x(2) = initial_velocity; % initial velocity % 2. Forward/Backward Euler time integration for n = 1 : nstep x = x + dt * ( A * x + b ); % Forward Euler y = (I - dt * A) \ (x + dt * b); % Backward Euler z = 0.5*(y + x); %Crank-Nicolson all_Euler_x (n,:) = x; % store solution all_Backward_y (n,:) = y; % store solution all_CN_z (n,:) = z; % store solution end % 3. Compute the exact solution % 3.1 parameters needed to compute the exact solution c_c = 2.0 * sqrt ( k * m ); zeta = c / c_c; omega_n = c_c / ( 2.0 * m ); c1 = ( initial_velocity + zeta * omega_n * initial_displacement ) ... / ( omega_n * sqrt ( 1.0 - zeta^2 ) ); c2 = initial_displacement; % 3.2 exact solution for n = 1 : nstep t = n * dt; all_exact_x (n,1) = exp( -zeta * omega_n * t ) ... * ( c1 * sin ( sqrt(1.0 - zeta^2) * omega_n * t ) ... + c2 * cos ( sqrt(1.0 - zeta^2) * omega_n * t ) ); end % 4. Plot results % 4.1 plot settings set(0 ,'DefaultTextInterpreter' , 'latex' ); set(0 ,'DefaultAxesUnits' , 'pixels'); myFontSize = 14; myFontName = 'Helvetica' ; figure('Position', [100 100 700 500]); hold on % 4.2 plot exact and approximate solution plot( all_t, all_exact_x(:,1), '--b', 'LineWidth', 2 ) plot( all_t, all_Euler_x(:,1), 'r' , 'LineWidth', 1 ) plot( all_t, all_Backward_y(:,1)); plot( all_t, all_CN_z(:,1)); % 4.3 more plot settings set(gca,'Box' , 'on' ,... 'Color' , 'white' ,... 'LineWidth' , 0.5 ,... 'TickDir' , 'in' ,... 'XMinorTick' , 'on' , ... 'YMinorTick' , 'on' , ... 'XGrid' , 'on' , ... 'YGrid' , 'on' , ... 'XMinorGrid' , 'off' , ... 'YMinorGrid' , 'off' , ... 'XLim' , [ 0 50] , ... 'YLim' , [ -1 1] , ... 'XTick' ,linspace(0,50,6),... 'YTick' ,linspace(-1,1,5),... 'TickLength' , [.02 .02] ,... 'FontName' , myFontName ,... 'FontSize' , myFontSize ); set(gcf, 'color', 'white'); xlabel('time (s)', 'FontName', myFontName, 'FontSize', myFontSize) ylabel('displacement (m)', 'FontName', myFontName, 'FontSize', myFontSize) hleg1 = legend('Exact','Approximation','Backward','C-N'); set(hleg1,'Location','SouthEast') 

Is there a reason as to why I'm not getting a plot out of my code? What can I do to fix this? Any and all help would be appreciated. I'd also like to add error bars if possible: how can I do this?