How to update an old mex code to work again with ode45?

Hello

Back in 94 the following mex code worked with ode45 (Matlab 4.2 on solaris 2.6). Although I have the feeling that this type of code has been completely deprecated several years ago, I wonder whether it is possible to translate it to a new version (matlab 2017a). (I know that a m-file would do the job but I have lots of similar codes ..).

/* Lorenz Equation (84)  - Simulation */

#include "mex.h"
/* #include "mymex.h" */

/* Define the system sizes here: */
#define NSTATES 12
#define NOUTPUTS 12
#define NCOEFFS 16

static Matrix *Coeffs[NCOEFFS];

#define a mxGetPr(Coeffs[0])
#define b mxGetPr(Coeffs[1])
#define F mxGetPr(Coeffs[2])
#define G mxGetPr(Coeffs[3])
#define X0 mxGetPr(Coeffs[4])
#define X1 mxGetPr(Coeffs[5])
#define X2 mxGetPr(Coeffs[6])
#define X01 mxGetPr(Coeffs[7])
#define X11 mxGetPr(Coeffs[8])
#define X21 mxGetPr(Coeffs[9])
#define X02 mxGetPr(Coeffs[10])
#define X12 mxGetPr(Coeffs[11])
#define X22 mxGetPr(Coeffs[12])
#define X03 mxGetPr(Coeffs[13])
#define X13 mxGetPr(Coeffs[14])
#define X23 mxGetPr(Coeffs[15])


static
#ifdef __STDC__
void init_conditions(double *x0)
#else
void init_conditions(x0)
double *x0;
#endif /* __STDC__ */

{
x0[0] = X0[0];
x0[1] = X1[0];
x0[2] = X2[0];
x0[3]  = X01[0];
x0[4]  = X11[0];
x0[5]  = X21[0];
x0[6]  = X02[0];
x0[7]  = X12[0];
x0[8]  = X22[0];
x0[9]  = X03[0];
x0[10] = X13[0];
x0[11] = X23[0];

}

static
#ifdef __STDC__
void derivatives(double t, double *x, double *u, double *dx)
#else
void derivatives(t,x,u,dx)
double t, *x, *u;       /* Input variables */
double *dx;             /* Output variable*/
#endif /* __STDC__ */

{

int i;

dx[0] = -x[1] * x[1] - x[2] * x[2] - a[0] * (x[0] - F[0]);
dx[1] = x[0] * x[1] - b[0] * x[0] * x[2] - x[1] - G[0];
dx[2] = b[0] * x[0] * x[1] + x[0] * x[2] - x[2];

for (i=0;i<=2;i++)
        {
        dx[3+i] = (-a[0]) * x[3+i] + (-2.0 * x[1]) * x[6+i] + (-2.0 * x[2]) * x[9+i];
        dx[6+i] = (x[1] - b[0] * x[2]) * x[3+i] + (x[0] - 1.0) * x[6+i] + (-b[0] * x[0]) * x[9+i];
        dx[9+i] = (b[0] * x[1] + x[2]) * x[3+i] + (b[0] * x[0]) * x[6+i] + (x[0] - 1.0) * x[9+i];
        }

}

static
#ifdef __STDC__
void outputs(double t, double *x, double *u, double *y)
#else
void outputs(t,x,u,y)
double t, *x, *u;       /* Input variables */
double *y;              /* Output variable*/
#endif /* __STDC__ */

{
int i;
for (i=0;i<NOUTPUTS;y[i]=x[i],i++);
}

/* #include "simulink.h" */

The function was called by issuing

[t,x,y]=ode45('lor84',tf,[],options,[],a,b,F,G,X0,X1,X2,X3,X4,X5,X6,X7,X8,X9,X10,X11);

Many thanks.

Ed

[James Tursa]

This code looks incomplete. You don't show Coeffs being filled with any values, and you don't show any code calling these functions. Absent any of that information, the only thing I can suggest is possibly changing this

static Matrix *Coeffs[NCOEFFS];

to this

static mxArray *Coeffs[NCOEFFS];

But even that is a guess since you don't show what Matrix is.

[Ed Mendes]

Many thanks. As I said the code worked back in 94 without a problem. Matrix was defined in the old mex.h that came with Matlab 4.2 and it is likely that mxArray is its successor. The c code, lor84.c, is called using ode45 as show in my first post. I think the big problem is simulink.h that does not give the support for such a function as it did back in 94.

[Ed Mendes]

I think I found the problem - "error: Pre-Simulink 2.0 S-functions no are longer supported. We suggest using one of the S-function templates and adding your algorithm code to it. #error Pre-Simulink 2.0 S-functions no are longer supported. We suggest using one of the S-function templates and adding your algorithm code to it."

It seems that I have to use a S-function c-template.

There is an example, csfunc.c, that can be useful. If you have an idea how to change it to simulate the system in my c-code above, please let me know. Many thanks.

[Ed Mendes]

There is another interesting function, vdpmex.c. The only thing that I could not figure out is how to pass the parameters as I did in the c-code above, that is, a,b,F,G and initial conditions. And of course, how to call it from a command line using ode45. many thanks.