//for call in matlab [x,p]=eilertme(tau,nstep);
#include"mex.h"
#include<stdlib.h>
#include"vector.h"
#include<math.h>
template<class A> void eiler(A & y, double tau, int nstep, A(*f)(A, double));

template<class A> void eiler(A & y, double tau, int nstep, A(*f)(A, double)) {
    for (int i = 0; i < nstep; i++)y = y + tau * f(y, i * tau);
}
Vector right(Vector y, double t);

Vector right(Vector y, double t) {
    Vector tmp(2);
    tmp(1) = y(2);
    tmp(2) = sin(y(1));
    return tmp;
}

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
    Vector y(2);
    double tau;
    double y1[100], y2[100];
    double *d1, *d2;
    int nstep;
    tau = *mxGetPr(prhs[0]);
    nstep = (int) * mxGetPr(prhs[1]);
    y(1) = 0.1;
    y(2) = 0.1;
    for (int i = 0; i < 100; i++) {
        eiler(y, tau, nstep, right);
        y1[i] = y(1);
        y2[i] = y(2);
    }
    plhs[0] = mxCreateDoubleMatrix(100, 1, mxREAL);
    plhs[1] = mxCreateDoubleMatrix(100, 1, mxREAL);

    d1 = mxGetPr(plhs[0]);
    d2 = mxGetPr(plhs[1]);
    for (i = 0; i < 100; i++) {
        d1[i] = y1[i];
        d2[i] = y2[i];
    }

}