function [nt,t,x,dxdp] = solve_map(x0,p,tfinal,xdim,kdim,solution_only)
     
% Initialize arrays
% The final xdim sensitivities are wrt the initial conditions

x=zeros(xdim,tfinal+1);
t=[0:tfinal];
nt=tfinal+1;
dxdp=zeros(xdim,kdim+xdim,tfinal+1);

if(solution_only == 0)
    for i=1:xdim
         dxdp(i,kdim+i,1)=1;
    end
end

x(:,1)=x0;
for it = 1:tfinal 
    x(:,it+1)=gvec(x(:,it),p);
end

if(solution_only == 0)
    for it = 1:tfinal 
        dxdp(:,:,it+1)=drhs_dparam(x(:,it),dxdp(:,:,it),p);
    end
end

end
