#include <math.h>
#include <stdio.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv.h>


int solow_de(double t, const double y[], double f[], void *params)
{
    double *p = (double *) params;
    double alpha = p[0];
    double a = p[1];
    double b = p[2];
    double k = y[0];

    /* dk/dt = a*k^alpha - b*k */
    f[0] = a*pow(k,alpha) - b*k;
    return GSL_SUCCESS;
}

int jac(double t, const double y[], double *dfdy, double dfdt[], void *params)
{
    double *p = (double *)params;
    double alpha = p[0];
    double a = p[1];
    double b = p[2];
    double k = y[0];
    
    gsl_matrix_view dfdy_mat 
        = gsl_matrix_view_array(dfdy, 1, 1);
    gsl_matrix * m = &dfdy_mat.matrix;

    double der = a*alpha*pow(k, alpha-1.0) - b;

    gsl_matrix_set (m, 0, 0, der);
    dfdt[0] = 0.0;
    return GSL_SUCCESS;
}

int main (void)
{
    double alpha = 1.0/3.0;
    double a = 1.0;
    double b = 0.5;
    double params[3];
    const gsl_odeiv_step_type * T = gsl_odeiv_step_rk8pd;

    gsl_odeiv_step * s  = gsl_odeiv_step_alloc (T, 1);
    gsl_odeiv_control * c  = gsl_odeiv_control_y_new (1e-9,1e-7);
    gsl_odeiv_evolve * e  = gsl_odeiv_evolve_alloc (1);

    params[0] = alpha;
    params[1] = a;
    params[2] = b;

    gsl_odeiv_system sys = {solow_de, jac, 1, params};

    double t = 0.0;
    double t1 = 30.0;
    double h = 1e-6;
    double y[1] = { 0.05 };

    while (t < t1)
    {
        int status = gsl_odeiv_evolve_apply (e, c, s,
                                        &sys, 
                                        &t, t1,
                                        &h, y);
        if (status != GSL_SUCCESS)
            break;
        printf ("%.5e %.5e\n", t, y[0]);
    }

    gsl_odeiv_evolve_free (e);
    gsl_odeiv_control_free (c);
    gsl_odeiv_step_free (s);
    return 0;
}