function [pvec,rho] = run_painsim_pde(patient_info, sim_info, drug_info, dose_times, pain_report_real)
% rate of decrease of pain without drugs or disease
relax_rate = patient_info(1);
% baseline pain, equilibrium for patient in absence of drugs
base_pain = patient_info(2);
% initial pain
init_pain = patient_info(3);
% noise strength
eps = patient_info(4);

% initial time
t0 = sim_info(1);
% end time
tmax = sim_info(2);
% time resolution
dt = sim_info(3);
% time vector
tvec = t0:dt:tmax;

% read in "real" patient data
report_times_real = pain_report_real(1,:);

% run pain simulation with one drug intervention
[tvec,pvec,rho] = pain_sim_drug1_pde(tvec,init_pain,relax_rate,base_pain,eps,drug_info,dose_times);

% report distribution of pain level at each time requested
[T_mesh,P_mesh] = meshgrid(tvec,pvec');
[T_interp,P_interp] = meshgrid(report_times_real,pvec');
rho = interp2(T_mesh,P_mesh,rho',T_interp,P_interp);


end