function [nt,t,x,dxdp] = solve_ode(x0,p,tfinal,xdim,kdim,solution_only)

    if xdim ~= length(x0)
        fprintf('solve_ode: Houston, we have a problem \n')
    end
    if kdim ~= length(p)
        fprintf('solve_ode: Houston, we have a problem \n')
    end

    % Set initial conditions for solution and for sensitivity computations
    % If solution_only = 1, then want only x0 
    if solution_only == 0
        xx0=[x0 zeros(xdim,kdim) eye(xdim)];
    end
    
    % Matlab routine ode45 solves vector systems of odes not matrix-valued odes
    % Reshape the sensitivity computations to a vector systems by unpacking column by column
    % to create a vector of length xdim*(1+kdim+xdim)
    if solution_only == 0
        y0=reshape(xx0,xdim*(1+kdim+xdim),1);
        options = odeset('RelTol',1E-6,'AbsTol', 1E-06*ones(xdim*(1+kdim+xdim),1));
    else
        y0 = x0;
        options = odeset('RelTol',1E-6,'AbsTol', 1E-06*ones(xdim,1));
    end
    
    % Solve system of odes,   
    %      xdim nonlinear odes
    %      kdim*xdim odes for sensitivity to parameters
    %      xdim*xdim odes for sensitivity to initial conditions

    % ode45 expects a handle to a function with two parameters t and y to create the rhs
    % We therefore use a "nested function" gsystem, which expects five parameters,
    %       t, y, p, xdim and kdim
    %  to create the rhs 
    %  Here y is the column vector of length xdim*(1+kdim+xdim) containing the solution 
    %  and sensitivity information, i.e.,
    %       y = [x; dxdp(:,1); dxdp(:,2); ... dxdp(:,kdim); ... dxdp(:,kdim+xdim)];

    [t,y]=ode45(@(t,y) gsystem(t,y,p,xdim,kdim,solution_only), [0 tfinal], y0, options);

    % Reverse row and column order of y as returned by ode45
    t=t';
    y=y';
    fprintf('size(t) = (%5i, %5i) \n', size(t))
    fprintf('size(y) = (%5i, %5i) \n', size(y))
    [ny,nt]=size(y);

    % Map solution and sensitivities for printing and plotting
    x(1:xdim,1:nt)=y(1:xdim,1:nt);
    
    dxdp = zeros(xdim,kdim+xdim,nt);
    if(solution_only == 0)
        for k=1:kdim+xdim
            dxdp(1:xdim,k,1:nt)=y(k*xdim+1:(k+1)*xdim,1:nt);
        end
    end
    
 end

