Runge-Kutta 2nd order of differential equations Could I improve this somehow (not numerically, just C-wise)?  Do I need to allocate the array since I am returning it and not just using it inside the function?

#include &lt;stdio.h&gt;
#include &lt;stdlib.h&gt;
#include &lt;math.h&gt;

/* Approximates a solution to a differential equation on the form: 
   y'(t) + ay(t) = x(t)
   y(0) = b 
*/
double* runge_kutta_2nd_order(double stepSize, double a, double b, double (*x) (double), double upto)
{
    int resultSize = ((int) (upto / stepSize)) + 2;
    double yt = b;
    double time;
    double k1,k2,ystar1,ystar2;
    int index = 1;

    //double *results = (double*) malloc(resultSize * (sizeof(double)));
    double *results = malloc(resultSize * sizeof(*results));
    if(results == NULL)
    {
        printf("\nCould not allocate memory. Exiting program.");
        exit(0);
    }

    results[0] = b;

    for(time = 0; time &lt; upto; time += stepSize) //&lt;=
    {   
        k1 = x(time) - a * yt;
        ystar1 = yt + stepSize * k1;
        k2 = x(time + stepSize) - a * ystar1;
        ystar2 = yt + (k1 + k2) / 2 * stepSize;
        yt = ystar2;
        results[index] = ystar2;
        index++;
    }
    return results;
}

void free_results(double **r)
{
    free(*r);
    *r = NULL;
}


double insignal(double t)
{
    return exp(t/2)*(sin(5*t) - 10*cos(5*t));
}

int main(void)
{
    int i;
    double *res = runge_kutta_2nd_order(0.01,-1,0,&amp;insignal,10);


    printf("\nRunge Kutta 2nd order approximation of the differential equation:");
    printf("\ny'(t) - y(t) = e^(t/2) * (sin(5t) - 10cos(5t))");
    printf("\ny(0) = 0");
    printf("\n0 &lt;= t &lt;= 10");

    for(i=0; i&lt;1001; i+=100){
        printf("\ni = %lf =&gt; y = ", 0.01*i);
        printf("%lf", res[i]);
    }
    printf("\n");

    free_results(&amp;res);


    return 0;
}
It is more usual to take a pointer to an already-allocated array, rather than allocate an array within the function itself.  This allows the calling function to choose the best kind of allocation.

When doing this you'd likely also take the number of indexes to fill as an input, deriving the step size from that:

/* resultsize must be &gt;= 2 */
void runge_kutta_2nd_order(double results[], size_t resultsize, double a, double b, double (*x) (double), double upto)
{
    const double stepSize = upto / (resultsize - 1);
    double yt = b;
    size_t index;

    results[0] = b;

    for (index = 1; index &lt; resultsize; index++)
    {
        const double time = index * stepSize;
        const double k1 = x(time) - a * yt;
        const double ystar1 = yt + stepSize * k1;
        const double k2 = x(time + stepSize) - a * ystar1;
        const double ystar2 = yt + (k1 + k2) / 2 * stepSize;
        yt = ystar2;
        results[index] = ystar2;
    }
}


So, for example, in this case you could use automatic allocation in main():

int main(void)
{
    int i;
    double res[1001];

    runge_kutta_2nd_order(res, 1001, -1, 0, &amp;insignal, 10);

    printf("\nRunge Kutta 2nd order approximation of the differential equation:");
    printf("\ny'(t) - y(t) = e^(t/2) * (sin(5t) - 10cos(5t))");
    printf("\ny(0) = 0");
    printf("\n0 &lt;= t &lt;= 10");

    for(i=0; i&lt;1001; i+=100){
        printf("\ni = %f =&gt; y = ", 0.01*i);
        printf("%f", res[i]);
    }
    printf("\n");

    return 0;
}