Metodo di Newton-Raphson multidimensionale per ODE

Buonasera, sto cercando di implementare il metodo di Newton Raphson nel caso multidimensionale così da risolvere ODE del secondo ordine. Ho trovato questo link che ho provato a seguire https://skill-lync.com/student-projects/Multivariate-Newton-Raphson-Solver-for-ODE-s-87655 ma non credo di averlo applicato nel modo giusto.
In particolare quando utilizzo la funzione che svolge il metodo di Newton questo non trova mai la soluzione. Dove sto sbagliando?
Nelle funzioni dy1f1 e simili sono contenute le derivate parziali per la Jacobiana. Le funzioni e f2 rappresentano le ODE di seconda specie, in particolare sarebbero y1' e y2'.
L'output del programma deve essere il file euleroimplicito.txt con 3 colonne: t, y1, y2 .

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <cmath>
#define SIZE 10
#define _USE_MATH_DEFINES
#define delta 0.0005
FILE *file;

typedef float vector [SIZE];

void euleroimplicito2(double, double, double, int, double, int);
int funzioni(int);
void condizioni(double*, double*);
void tempi(double*, double*);
double f1(double, double, double, int);
double f2(double, double, double, int);
double dy1f1(double, double, int, double);
double dy2f1(double, double, int, double);
double dy1f2(double, double, int, double, double);
double dy2f2(double, double, int, double, double);
double newton2(double, double, double, double, int, int, vector);

int main() {
	double t0, tf, x0, v0, h;
	int f,N, i;
	f=funzioni(i);
	tempi(&t0, &tf);
    condizioni(&x0, &v0);
    printf("inserisci il passo h:  ");
	scanf("%lf", &h);
	N=(int)(tf-t0)/h+1;
    euleroimplicito2(t0,x0,v0,f,h,N);
}

void euleroimplicito2(double t0, double x0, double v0, int f, double h, int N){
	double t=t0;
	vector u1, u2, b;
	int j;
	FILE* file;
	file=fopen("euleroimplicito.txt","w");
	u1[0]=x0;
	u2[0]=v0;
	fprintf(file, "\n%lf %lf %lf",t,u1[0],u2[0]);
	for(j=1; j<=N; j++){
		newton2(u1[j-1],u2[j-1],h,t,f,j, b);
		fprintf(file,"\n%lf %lf %lf",t, b[0],b[1]);
		u1[j]=b[0];
		u2[j]=b[1];
		t=t0+h*j;
		}
	fclose(file);
}

double newton2(double x0, double v0, double h, double t, int f, int j, vector b)
{
	int k=1, q=0, N=100;
	vector rapp;
	double J[2][2], det;
	J[0][0]=dy1f1(x0,v0,f,h);
	J[0][1]=dy2f1(x0,v0,f,h);
	J[1][0]=dy1f2(x0,v0,f,t,h);
	J[1][1]=dy2f2(x0,v0,f,t,h);
	det=(J[0][0]*J[1][1])-(J[0][1]*J[1][0]);
	b[0]=0; b[1]=0;
	if(det==0) {printf("La matrice jacobiana associata all'ODE non è invertibile, non è possibile continuare"); exit(0);}
	double w1, w2, err=1e-3;
	do{
	    	J[0][0]=dy1f1(x0,v0,f,h);
         	J[0][1]=dy2f1(x0,v0,f,h);
        	J[1][0]=dy1f2(x0,v0,f,t,h);
         	J[1][1]=dy2f2(x0,v0,f,t,h);
         	det=(J[0][0]*J[1][1])-(J[0][1]*J[1][0]);
        	if(det==0) {printf("La matrice jacobiana associata all'ODE non è invertibile, non è possibile continuare"); exit(0);}
		//rapp[0]=J11/det[f1(x0,v0)]-J01/det*[f2(x0,v0)]
		rapp[0]=((J[1][1]/det)*f1(t,x0,v0,f))-((J[0][1]/det)*f2(t,x0,v0,f));
		//rapp[1]=-J10/det*[f1(x0,v0)]+J00/det*[f2(x0,v0)]
		rapp[1]=((-J[1][0]/det)*f1(t,x0,v0,f))+((J[0][0]/det)*f2(t,x0,v0,f));
		w1=x0-rapp[0];
		w2=v0-rapp[1];
		if((fabs(rapp[0])<err) && q==0)
		{
			b[0]=w1;
			b[1]=w2;
			q++;
		}
		x0=w1;
		v0=w2;
		k++;
	}while(k<=N && q==0);
	if (q==0) printf(" The required solution for EB does not converge or iterations are insufficient\n");
}

int funzioni(int i){
		int f;
		printf("\nLe equazioni che si possono risolvere sono: \n1)Oscillatore armonico \n2)Oscillatore armonico smorzato \n3)Oscillatore armonico forzato \n4)Pendolo con attrito \n5)Modello Lotka-Volterra \n6)Van der Pol \n");
	    do{ printf("\nScegliere il numero dell'equazione che si vuole risolvere:  ");
		scanf("%d", &f);
	    } while (f<1 || f>6);
    return f;
}

void condizioni(double* x0, double* v0){
	printf("Inserisci x0:  ");
	scanf("%lf", x0);
	printf("Inserisci v0:  ");
	scanf("%lf", v0);
}

void tempi(double* t0, double* tf){
	printf("Inserisci t0:  ");
	scanf("%lf", t0);
	printf("Inserisci tf:  ");
	scanf("%lf", tf);
}

double f1(double t, double y1, double y2, int f){  //y1=y, y2=y'
	switch(f){
		case 1:
		case 2:
		case 3:
		case 4:
			return y2;
		case 5:
			return 0.4*y1-0.4*y1*y2; //prede //x0=100
			break;
		case 6:
			return y2-y1*y1*y1+y1;
			break;
	}
}

double dy2f1(double y1, double y2, int f, double h){
	switch(f){
		case 1:
		case 2:
		case 3:
		case 4:
			return -h;
		case 5:
		    return h*0.4*y1;
		case 6:
			return -h;
    }
}

double dy1f1(double y1, double y2, int f, double h){
	switch(f){
		case 1:
		case 2:
		case 3:
		case 4:
			return 1;
		case 5:
			return 1-0.4*h+0.4*y2*h;
		case 6:
			return 1+3*h*y1*y1-h;
    }
}

double f2(double t, double y1, double y2, int f){
	switch(f){
		case 1:
			return -y1;
			break;
		case 2:
			return -9*y1-2*y2;
			break;
		case 3:
			return -4*y1+t*t;
			break;
		case 4:
			return -y2-sin(y1);
			break;
		case 5:
			return -2*y2+0.09*y1*y2; //predatori //y0=8 h=0.01
			break;
		case 6:
			return -y1;
			break;
	}
}

double dy2f2(double y1, double y2, int f, double t, double h){
	switch(f){
		case 1:
		return 1;
		case 2:
			return 1+2*h;
		case 3:
			return 1;
		case 4:
			return 1+h;
		case 5:
			return 1+2*h-0.09*h*y1;
		case 6:
		    return 1;
	}
}

double dy1f2(double y1, double y2, int f, double t, double h){
	switch(f){
		case 1:
		  return h;
		case 2:
		  return 9*h;
		case 3:
			return 4*h;
		case 4:
			return h*cos(y1);
		case 5:
		    return -h*0.09*y2;
		case 6:
			return h;
	}
}