Plotting multiple functions with changing constants

I would like to plot the equation for rho as a function of t on two graphs, one with beta=50 and one with beta=60. On each graph, I want to be able to plot each of the 5 curves created by the changing constant v0. Also, I would like for each curve to only be in the range provided. How would I be able to to this?

[Image Analyst]

Try this

% Initialization steps.
clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 15;
fprintf('Beginning to run %s.m ...\n', mfilename);

g = 9.8;
legendStrings = cell(5, 1);
loopCounter = 1;
for v0 = 30 : 10 : 70
    numElements = 1980; % HDTV resolution
    p1 = subplot(2, 1, 1);
    beta = 50;
    % Construct t
    t = linspace(0, 2 * v0 * sind(beta)/g, numElements);
    % Get rho
    rho = ComputeRho(beta, v0, t);
    plot(t, rho, '-', 'LineWidth', 2);
    grid on;
    hold on;
    legendStrings{loopCounter} = sprintf('v0 = %d', v0);
    title('Beta = 50')
    xlabel('t');
    ylabel('rho')
    drawnow;

    p2 = subplot(2, 1, 2);
    beta = 60;
    % Construct t
    t = linspace(0, 2 * v0 * sind(beta)/g, numElements);
    % Get rho
    rho = ComputeRho(beta, v0, t);
    plot(t, rho, '-', 'LineWidth', 2);
    grid on;
    hold on;
    title('Beta = 60')
    xlabel('t');
    ylabel('rho')
    loopCounter = loopCounter + 1;
    drawnow;
end
legend(p1, legendStrings, 'Location', 'north');
legend(p2, legendStrings, 'Location', 'north');

function rho = ComputeRho(beta, v0, t)
g = 9.8;
term1 = (v0^2 * cosd(beta).^2) / g;
term2 = -g*t ./ (v0  * cosd(beta));
term3 = sqrt((1 + (term2 + tand(beta)).^2).^3);
rho = term1 .* term3;
end

[Eli Wolkenstein]

Thank you