(Finally!) I figured out a solution using function handles: The following works when solving a scalar pde problem. See below for the pde system case... Scalar pde: 1) Create a function for C (or A, or F, etc) which uses the additional parameters C = C_fun(p,t,u,time,additional_parameters)
2) In the code calling assempde, create function handles for the coefficients eg: C_h = @(p,t,u,time)C_fun(p,t,u,time,additional_params); A_h = @(p,t,u,time)A_fun(p,t,u,time,additional_params); bc_h = @(p,t,u,time)bc_fun(p,t,u,time,additional_params); f_h = @(p,t,u,time)f_fun(p,t,u,time,additional_params);
3) Call assempde using the above function handles: u = assempde(bc_h,p,e,t,C_h,A_h,f_h)
System of pde's: If the above steps are used when you have a system of PDE's, then an error is given by assema.m saying that the number of rows in C is incorrect. This is because the number of rows N is defined by the number of rows in f. However, there's a bug in this code: the number of rows in f is evaluated from the number of rows in input argument for f (line 100 in assema.m). Hence, if f_h in the above example contains a single row (as is the case for a function handle) then N is set to 1. The actual evaluation of the f in assema.m occurs at the end of the program (lines 234 and down), starting with f=pdetfxpd(p,t,uu,time,f); If this is section of the code is simply moved above line 100, then f has number of rows given to it by the evaluation of the input argument for f (which should be N) and the bug is fixed.