Matlab/ Octave: Gráfica con plot( ... )

Uno de mis primeros ejercicios sobre métodos numéricos, a continuación un código que abarca algunos puntos sobre la programación en Matlab.
El código debe ser ejecutado en Matlab o en Octave (solución alternativa gratuita a Matlab). La gráfica fue hecha con Octave.
Puedes copiar el archivo sin incluir las líneas y guardarlo con extensión .m

--------------------------------------------- INICIO DE ARCHIVO --------------------------------

% Autor: Juan Carlos Miranda
% juancarlosmiranda81@gmail.com.
% www.otroblogdetecnologias.blogspot.com
% 03/12/2015
%
% Ejercicio en Matlab
% -----------------------
% En este ejemplo, se utilizan varios aspectos relacionados a la
% programación en Matlab, en base a un problema, se muestran varios
% items útiles en la resolución de problemas numéricos.

% Los aspectos que se desarrollan son:
%
% -Definición de variables.
% -Definición de vectores.
% -Uso de ciclos.
% -Ploteado de gráficas con una variable y con dos variables.
%
%
% Planteamiento del problema
% ---------------------------
% Dada la ecuación diferencial del aceleramiento de una masa en caída libre por, dv/dt=g-(cv/m) donde v es la velocidad en m/s, g es la constante gravitacional que tiene un valor
% de 9.8 m/s 2 , m es la masa con un valor de 68.1 kg y c es el coeficiente de
% resistencia con una valor de 12.5 kg/s.

% En la solución, se expondrán dos gráficas superpuestas, una con el
% resultado de la ecuación analítica y otro con el resultado del método de
% aproximación utilizando la fórmula de Euler.

% Inicializacion de variables de entorno
clear all;
close all;
clc;

% Datos del cuerpo en caida libre
g=9.8; %gravedad
c=12.5; % coeficiente de rozamiento
m=68.1; %masa

%Rango de tiempo para el cáculo.
ti=0;
tiempoFinal=14;

%paso de avance, refleja la cantidad que se incrementará en la Ecuación de
%aproximación por el método de Euler.
paso=0.5;
vi=0;  % velocidad con la que arranca el cuerpo.
va=0;  % definición de la velocidad por ecuación analítica
vti=0; % velocidad por método de Euler.

%inicializacion de vectores variables
tiempo=zeros(0,1,'double');
velocidadi=zeros(0,1,'double');
velocidada=zeros(0,1,'double');

%indice para el guardado en vectores
index=1;

% Impresión de cabecera
str = fprintf('-------------------------------------------------------\n');
str = fprintf('Ti; Vi ; Ti+h; VEuler(Ti); VAnalitica \n');   
str = fprintf('-------------------------------------------------------\n');

% se procede a calcular por cada unidad de tiempo luego se guarda en
% vectores que serán utilizados para la gráfica.
while ti < tiempoFinal
    %almacenamiento en vector de velocidad Euler y tiempo
    velocidadi(index)=vti;
    tiempo(index)=ti;

    %-----Ecuación del metodo Euler
    vti = vi+(g-((c/m)*(vi)))*(ti+paso);
    % Almacenamiento para el cálculo de velocidad en t+paso
    vi=vti;       
    %----------------------------
   
    %----Ecuación Analitica
    va=(g*m/c)*(1-exp((-c/m)*ti));
    %---guarda en vector resultado de ecuación analitica
    velocidada(index)=va;
   

    % Impresion de tabla
    str = fprintf('%d) %2.2i seg| %3.2f km/seg| %3.2f km/seg| \n', index,ti,vti,va);

    % incrementa el paso del indice de guardado de vectores
    index=index+1;
    %----se incrementa el paso del tiempo
    ti=ti+paso;

end

% gráfica con los vectores previamente cargados
plot(tiempo,velocidadi,'r');
%anexo de una segunda gráfica dentro de la figura, de manera a superponer
%las funciones.
figure(1),hold on,plot(tiempo,velocidada,'g');
% Definición de valores para representación del gráfico
set(gca,'fontsize',22);
xlabel('t=tiempo');
ylabel('Velocidad');
legend('Sol. Analitica','Sol. Numerica');
title('Ecuación de velocidad','FontSize',12);


--------------------------------------------- FIN DE ARCHIVO --------------------------------