A determinação da taxa mensal de juros pela determinação da raiz real positiva da função

[JuliaAlvesGuirado]

Olá, preciso de ajuda com uma atividade de linguagem c. O programa precisa ler o valor de certo bem, o número de prestações no qual foi pago, e o valor das prestações. A partir destes valores precisa calcular a taxa de juros pela raiz da função f(j)=somatório-(valor do bem/valor da prestação). Somatório é 1-(1/(1+j))^n/1-(1/(1+j)) quando j>0, sendo j a taxa de juros e n numero de prestações. Para achar a raiz da função é necessário usa o método de muller. Porém, o juros da astronômico, e já mexi em tudo e não sai. Se alguém poder ajudar a achar o erro.

O código é:

(ignorem as outras funções potencia, valor absoluto, erro relativo e raiz quadrada, é que não pode usar o math.h. O importante esta em negrito)

#include <stdio.h>
#include <stdlib.h>
#define TOL 0.000000001

double potencia(double x, int n)
{
    int k;
    double pot=1;

    for(k = 1; k <= n; k++)
    {
        pot = pot*x;
    }

    return pot;
}

double valor_absoluto(double x)
{
    if(x<0)
    {
      x=-1*x;
    }

    return x;

}

double erro_relativo(double x, double y)
{
    double erro=1, valor;
    if (y!=0)
        {
        valor=(x-y)/y;
        erro=valor_absoluto (valor);
        }
    else
        {
        erro=1;
        }
return erro;
}

double raiz_quadrada(double x)
{
    double esp=1, x0, x1;

    x0=x;

    while(esp>=TOL)
    {
    x1=x0-(((x0*x0)-x)/(2*x0));
    esp=valor_absoluto(((x1-x0)/x1));
    x0=x1;
    }
    return x1;
}

double f(double j,double vbem, int nprest, double prest)
{
    double somatorio, funcao;

    if (j==0)
    {
        somatorio=nprest;
        funcao=somatorio-(vbem/prest);
    }
    else
    {
        somatorio=(1-(potencia(1/(1+j),nprest)))/(1-(1/(1+j)));
        funcao=somatorio-(vbem/prest);
    }

    return funcao;
}

double juros (double vbem, int nprest, double prest)
{
double x0=0, x1, fx0, fx1, x2, fx2, fx3, h1, h2, q2, a, b, c, x3, esp, fx=1, resp, qlinha, qlinhalinha;

x1=0.3*(prest/(vbem-prest));
x2=prest/(vbem-prest);
fx0=f(x0,vbem,nprest,prest);
fx1=f(x1,vbem,nprest,prest);
fx2=f(x2,vbem,nprest,prest);
h1=x1-x0;
h2=x2-x1;
q2=h2/h1;
a=(q2*fx2)-(q2*(1+q2)*fx1)+((potencia(q2,2))*fx0);
b=((2*q2+1)*fx2)-((potencia((1+q2),2))*fx1)-((potencia(q2,2))*fx0);
c=(1+q2)*fx2;
qlinha=(-2*c)/(b+(raiz_quadrada(b*b-(4*a*c))));
qlinhalinha=(-2*c)/(b-(raiz_quadrada(b*b-(4*a*c))));
if (valor_absoluto(qlinha)>valor_absoluto(qlinhalinha))
{
    x3=x2+(h1*qlinhalinha);
}
else
{
    x3=x2+(h1*qlinha);
}
fx3=f(x3,vbem,nprest,prest);
esp=erro_relativo(x3, x2);
x0=x1;
x1=x2;
x2=x3;
resp=x3;

if (fx0==0)
{
    fx=fx0;
    resp=x0;
}

if (fx1==0)
{
    fx=fx1;
    resp=x1;
}

if (fx2==0)
{
    fx=fx2;
    resp=x2;
}

while (esp>=TOL && fx!=0)
{
    fx0=f(x0,vbem,nprest,prest);
    fx1=f(x1,vbem,nprest,prest);
    fx2=f(x2,vbem,nprest,prest);
    h1=x1-x0;
    h2=x2-x1;
    q2=h2/h1;
    a=(q2*fx2)-(q2*(1+q2)*fx1)+((potencia(q2,2))*fx0);
    b=((2*q2+1)*fx2)-((potencia((1+q2),2))*fx1)-((potencia(q2,2))*fx0);
    c=(1+q2)*fx2;
    qlinha=(-2*c)/(b+(raiz_quadrada(b*b-(4*a*c))));
    qlinhalinha=(-2*c)/(b-(raiz_quadrada(b*b-(4*a*c))));
    if (valor_absoluto(qlinha)>valor_absoluto(qlinhalinha))
    {
        x3=x2+(h1*qlinhalinha);
    }
    else
    {
        x3=x2+(h1*qlinha);
    }

    fx3=f(x3,vbem,nprest,prest);
    esp=erro_relativo(x3, x2);
    x0=x1;
    x1=x2;
    x2=x3;
    resp=x3;

    if (fx0==0)
    {
        fx=fx0;
        resp=x0;
    }

    if (fx1==0)
    {
        fx=fx1;
        resp=x1;
    }

    if (fx2==0)
    {
        fx=fx2;
        resp=x2;
    }
}
return resp;
}

int main()
{
    int nprest=1;
    double vbem=1,prest=1;

    while (vbem>0)
    {
        printf("Digite valor do bem,numero de prestacoes e valor das prestacoes=");
        scanf("%lf %lf %lf", &vbem, &nprest, &prest);
        printf("juros= %lf\n", juros(vbem,nprest,prest));
    }

    system("PAUSE");
    return 0;
}