function runMatlabCellML1_4(display)
clear all;
clc;
tic;
global modelname Y constantNames constantValues computedNames computedValues stateNames

%%%Set up some common variables needing to be changed between models
time = 'time'; %Some models use 't', some use 'time'
lengthToRun = 0:1:12000; %Length of time to run the model for
suffix = '_COR_CuminEdit'; %Mark at end of saved file to show it's being used in this function

if nargin == 0
    display =1;
end
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Read in file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[filename path] = uigetfile('*.*');
t=textread([path,filename],'%s','delimiter','\r');

if(display)
    disp(['Using model: ',filename(1:strfind(filename,'.')-1),'     (time elapsed = ',num2str(toc),'s)']);
end

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Find variables & constants
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Set up to find where things are
textToFind = {'Initial conditions';'State variables';...
                'Constants';'Computed variables';'State variables'};
textPosition = zeros(2,length(textToFind));
%%%Find where things are

%%%Fill in blank lines so when saves, gap remains
for i=1:length(t)
    if strcmp(t{i},'')
        t{i} = '%';
    end
end

%%%Find start and end of sections
for i=1:length(t)
    for t2f = 1:length(textToFind);
        if (strfind(t{i}, textToFind{t2f}))
            textPosition(1,t2f) = i+3;
            for j=i+3:length(t)
                if( strfind(t{j},'%-') )
                    textPosition(2,t2f) = j-2;
                    break;
                end
            end
        end
    end
end

%%%INITIAL CONDITIONS
%%%Get the initial conditions to run the model with
Y = (t{textPosition(1,1)});
Y = eval(Y(6:length(Y)));
      
%%%STATE VARIABLES
%%%Store the constants for later use
if( (textPosition(1,5) > 0) && (textPosition(2,5)>0) )
    stateCount=1;   
    stateNames={};
    for i=textPosition(1,5):textPosition(2,5)
        lineText = t{i};
        %%Save constant and 'original' value for later if needed.
        stateNames{stateCount} = strtrim(lineText(strfind(lineText,':')+1:strfind(lineText,'(')-1));
        stateCount = stateCount+1;
    end
end

%%%CONSTANTS
%%%Store the constants for later use
if( (textPosition(1,3) > 0) && (textPosition(2,3)>0) )
    constantCount=1;
    constantNames={};
    constantValues=[]; %%Need to initialise the matrix

    for i=textPosition(1,3):textPosition(2,3)
        lineText = t{i};
        %%Save constant and 'original' value for later if needed.
        constantNames{constantCount} = strtrim(lineText(1:strfind(lineText,'=')-1));
        tempValue = str2num(lineText(strfind(lineText,'=')+1:...
                                                 strfind(lineText,';')-1));

        constantValues(constantCount) = tempValue;
        constantCount = constantCount+1;
    end
end

%%%COMPUTED VARIABLES
%%%Store computed variable names
if( (textPosition(1,4) > 0) && (textPosition(2,4)>0) )
    computedCount=1;
    computedNames={};
    for i=textPosition(1,4):textPosition(2,4)
       lineText = t{i};
       varName = lineText(2:strfind(lineText,'(')-2);
       computedNames{computedCount} = strtrim(varName);
       computedCount = computedCount+1;
    end
end

if(display)
    disp(['All varaibales read in.   (time elapsed = ',num2str(toc),'s)']);
end

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Edit the file to use later
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%COMPUTED VARS - so that we can easily plot them (or use later) if we want
%%%Working from the bottom up so that we don't 'cross over'
%%%Add variables to 'computed'
if( (textPosition(1,4) > 0) && (textPosition(2,4)>0) )
    for i=1:length(computedNames)
        lineText = strcat('computed(:,',num2str(i),')=',computedNames{i},';');
        t=[t;lineText];
    end
    %%% Initialise 'computed'
    t{14} = 'computed = [];';
    %%% add 'computed' to the function line
    t{13} = ['function [dY,computed] = ',filename(1:length(filename)-2),suffix,'(',time,',Y)'];
else
    %%% Don't initialise 'computed'
    t{14} = '%';
    %%% Don't add 'computed' to the function line
    t{13} = ['function [dY] = ',filename(1:length(filename)-2),suffix,'(',time,',Y)'];
end

%%%Write out new file
writeascii([path,filename(1:strfind(filename,'.')-1),suffix,'.m'],t,'%r\n');

%%Read in new file and use this one
filename=[filename(1:length(filename)-2),suffix,'.m'];
t=textread([path,filename],'%s','delimiter','\r');
    

modelname = filename(1:length(filename)-2);

if(display)
    disp(['File ready to be used   (time elapsed = ',num2str(toc),'s)']);
end


%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Compute the model and variables
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
time=lengthToRun;
[tData, YData] = runModel(modelname, time);

%%Save all to workspace
assignin('base','Y',Y);
for i=1:length(constantNames)
    assignin('base',constantNames{i},constantValues(i));
end
for i=1:length(computedNames)
    assignin('base',computedNames{i},computedValues(:,i));
end
for i=1:length(stateNames)
    assignin('base',stateNames{i},YData(:,i));
end
assignin('base','t',tData);

%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Delete the file once used
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
delete([path filename]);
toc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end %%%End of function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%Run the model
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [tData, YData] = runModel(modelname, time)
    global Y constantNames constantValues computedNames computedValues stateNames    

    fhandle = str2func(modelname);
    [tData, YData] = ode45(fhandle, time, Y);
    if(length(computedNames)>0)
        for i=1:length(tData)
            eval(['[dY,computedValues(i,:)] = ',modelname,'(tData(i),YData(i,:));']);
        end
        %%%Save computed variables as their names
        for i=1:length(computedNames)
            eval([computedNames{i},'=[',num2str(computedValues(:,i)'),'];']);
        end
        %%%Save computed variables as their names
        for i=1:length(stateNames)
            eval([stateNames{i},'=[',num2str(YData(:,i)'),'];']);
        end
    end
    %%%Save the constants as their names
    for i=1:length(constantValues)
        eval([constantNames{i},'=',num2str(constantValues(i)),';']);
    end

end