function mm_interface(gparam,gvec,qoi,cp,solution_only,NextGen,x0,p,imap,R0_only, DIR)
%
% Runs maple procedures from matlab
tic;
fprintf('\nmm_interface \n')
qoi_only=0;
if qoi_only==1
    fprintf('\n\nWARNING: SENSAI will produce QoI.m and cp.m ONLY \n\n')
end

xdim=gparam(1);
kdim=gparam(2);

fprintf('mm_interface: xdim =  %3i \n', xdim)
fprintf('mm_interface: kdim =  %3i \n', kdim)
for i=1:xdim
    %fprintf('mm_interface: g(%3i) is a string (true/false=1/0) %3i \n', i,ischar(g(i,:)))
    fprintf('mm_interface: g(%3i) =  %s \n', i,gvec(i,1:length(gvec(i,:))) )
end
fprintf('mm_interface: cp =  %s \n', cp)

gtotal=gvec(1,1:length(gvec(1,:)));
for i=2:xdim
    gtotal=strcat(gtotal,';',gvec(i,1:length(gvec(i,:))));
end
gtotal=strcat(' " ',gtotal,' " ')
%fprintf('mm_interface: gtotal is a string (true/false=1/0) %3i \n', ischar(gtotal))

fprintf('mm_interface: call dummy \n')
maple('read "MapleRoutines/dummy.mpl" ');
maple('dummy()');

% gvec maple files
fprintf('mm_interface: call gvec_top \n')
maple('read "MapleRoutines/gvec_top.mpl" ');
maple('gvec_top()');

fprintf('mm_interface: call gvec_middle \n')
maple('read "MapleRoutines/gvec_middle.mpl" ');
maple('gvec_middle', gparam,gtotal);

fprintf('mm_interface: call gvec_bottom \n')
maple('read "MapleRoutines/gvec_bottom.mpl" ');
maple('gvec_bottom', gparam);

if NextGen(1) ~= 0    
    
    % Create Next Generation Matrix using Maple
   fprintf('mm_interface: call r0_top \n')
   maple('read "MapleRoutines/r0_top.mpl" ');
   maple('r0_top()');

    if R0_only == 1      
       % Sensitivities of R0 are too difficult for SENSAI, compute R0 numerically
       fprintf('mm_interface: call r0_only_middle \n')
       maple('read "MapleRoutines/r0_only_middle.mpl" ');
       maple('r0_only_middle',NextGen',gparam,gtotal,x0,p,imap);
    else
       % Compute symbolic R0 and sensitivities
       fprintf('mm_interface: call r0_middle \n')
       maple('read "MapleRoutines/r0_middle.mpl" ');
       maple('r0_middle',NextGen',gparam,gtotal,x0,p,imap);
    end
    R0_time = toc;
    if R0_time > 60 
        fprintf(['R0 took ',num2str(R0_time/60),' MINUTES to compute.\n']);
    else
        fprintf(['R0 took ',num2str(R0_time),' SECONDS to compute.\n']);
    end
   fprintf('mm_interface: call r0_bottom \n')
   maple('read "MapleRoutines/r0_bottom.mpl" ');
   maple('r0_bottom',NextGen',gparam);

end
    
if(solution_only == 0)
    
    if(qoi_only==0)      
        
        % dgvec_dxvec maple files  
        fprintf('mm_interface: call dgvec_dxvec_top \n')
        maple('read "MapleRoutines/dgvec_dxvec_top.mpl" ');
        maple('dgvec_dxvec_top()');
        
        fprintf('mm_interface: call dgvec_dxvec_middle \n')
        maple('read "MapleRoutines/dgvec_dxvec_middle.mpl" ');
        maple('dgvec_dxvec_middle', gparam,gtotal);
        
        fprintf('mm_interface: call dgvec_dxvec_bottom \n')
        maple('read "MapleRoutines/dgvec_dxvec_bottom.mpl" ');
        maple('dgvec_dxvec_bottom', gparam);

        % dgvec_dparam maple files
        fprintf('mm_interface: call dgvec_dparam_top \n')
        maple('read "MapleRoutines/dgvec_dparam_top.mpl" ');
        maple('dgvec_dparam_top()');
        
        fprintf('mm_interface: call dgvec_dparam_middle \n')
        maple('read "MapleRoutines/dgvec_dparam_middle.mpl" ');
        maple('dgvec_dparam_middle', gparam,gtotal);
        
        fprintf('mm_interface: call dgvec_dparam_bottom \n')
        maple('read "MapleRoutines/dgvec_dparam_bottom.mpl" ');
        maple('dgvec_dparam_bottom', gparam);
        
    end
    
    % qoi maple files --- caution: qoi takes the derivative in it- make 2
    % maple files if want to keep qoi!
    fprintf('mm_interface: qoi_top \n')
    maple('read "MapleRoutines/qoi_top.mpl" ');
    maple('qoi_top()');
    
    fprintf('mm_interface: call qoi_middle \n')
    maple('read "MapleRoutines/qoi_middle.mpl" ');
    maple('qoi_middle', gparam,qoi);
    
    fprintf('mm_interface: call qoi_bottom \n')
    maple('read "MapleRoutines/qoi_bottom.mpl" ');
    maple('qoi_bottom',gparam);
    
    % cp maple files
    fprintf('mm_interface: cp_top \n')
    maple('read "MapleRoutines/cp_top.mpl" ');
    maple('cp_top()');
    
    fprintf('mm_interface: call cp_middle \n')
    maple('read "MapleRoutines/cp_middle.mpl" ');
    maple('cp_middle', gparam,cp);
    
    fprintf('mm_interface: call cp_bottom \n')
    maple('read "MapleRoutines/cp_bottom.mpl" ');
    maple('cp_bottom',gparam);
    
%         % dcp_dparam maple files
%         fprintf('mm_interface: call dcp_dparam_top \n')
%         maple('read "MapleRoutines/dcp_dparam_top.mpl" ');
%         maple('dcp_dparam_top()');
%     
%         fprintf('mm_interface: call dcp_dparam_middle \n')
%         maple('read "MapleRoutines/dcp_dparam_middle.mpl" ');
%         maple('dcp_dparam_middle', gparam,cp);
%     
%         fprintf('mm_interface: call dcp_dparam_bottom \n')
%         maple('read "MapleRoutines/dcp_dparam_bottom.mpl" ');
%         maple('dcp_dparam_bottom',gparam);
    
end

toc;
end
