Maxima is mostly an algebraic system , however using various simplifications numerical calculations are also possible ; even those that have a tendency to cause numerical overflow. Yes numerical overflow is possible in maxima unless various precautions are taken. This is the maxima equivalent to the fractint formula d2jbMandelbrot that I posted recently . You may not need to compile() the function , I tend to do this as I'm always assuming that a speed increase will result. The file iterdiff2c.wxm /* [wxMaxima batch file version 1] [ DO NOT EDIT BY HAND! ]*/ /* [ Created with wxMaxima version 11.08.0 ] */ /* [wxMaxima: comment start ] . iterdiff2c.wxm (c) Copyright 2017 , sciwise@ihug.co.nz Methodology for expressing differential equations in terms of an iterated function. In particular a first order differential equation and the derivative thereof . Note that f(z,c) might be any iterated function . . [wxMaxima: comment end ] */ /* [wxMaxima: comment start ] ' Differential equation we're examining . diff(f(z),z) + 3*f(z)*z = sin(z) for an iterated function z = f(z) , g(z) =f(f(z)) diff(g(z),z) + 3*g(z)*f(z) = sin(f(z)) . [wxMaxima: comment end ] */ /* [wxMaxima: input start ] */ kill(all)$ eq1:diff(g(z,c),z) + 3*g(z,c)*f(z,c) - sin(f(z,c))$ expand(eq1); expand(diff(eq1,z)) ; /* [wxMaxima: input end ] */ /* [wxMaxima: comment start ] For the iterated Mandelbrot fractal function . [wxMaxima: comment end ] */ /* [wxMaxima: input start ] */ f(z,c):=z*z+c$ g(z,c):=f(f(z,c),c)$ f(z,c); g(z,c); /* [wxMaxima: input end ] */ /* [wxMaxima: comment start ] [wxMaxima: comment end ] */ /* [wxMaxima: input start ] */ eq1:diff(g(z,c),z) + 3*g(z,c)*f(z,c) - sin(f(z,c))$ ed:expand(eq1)$ edp:expand(diff(eq1,z))$ /* [wxMaxima: input end ] */ /* [wxMaxima: comment start ] [wxMaxima: comment end ] */ /* [wxMaxima: input start ] */ ed; edp; /* [wxMaxima: input end ] */ /* [wxMaxima: comment start ] This may well simplify the construction of fractint formulas for this type of equation . [wxMaxima: comment end ] */ /* [wxMaxima: comment start ] z=Pixel c=Pixel : ed = -sin(z^2+c)+3*z^6+9*c*z^4+4*z^3+9*c^2*z^2+3*c*z^2+4*c*z+3*c^3+3*c^2 edp = -2*z*cos(z^2+c)+18*z^5+36*c*z^3+12*z^2+18*c^2*z+6*c*z+4*c [wxMaxima: comment end ] */ /* [wxMaxima: input start ] */ kill(all)$ /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ ratprint : false$ /* [wxMaxima: input end ] */ /* [wxMaxima: comment start ] Iterated function , Differential and Newton's method. [wxMaxima: comment end ] */ /* [wxMaxima: input start ] */ iterdiffn(c,n):=block([z,ed,edp,mag,j,i], mag:0.001, j:0, z:c, ed:1, for i:0 while( i < n) and ( abs(ed) > mag ) do( z:expand(z), z:radcan(z), z:realpart(z) + %i*imagpart(z), z:ev(z,numer), ed : -sin(z^2+c)+3*z^6+9*c*z^4+4*z^3+9*c^2*z^2+3*c*z^2+4*c*z+3*c^3+3*c^2, ed:expand(ed), ed:radcan(ed), ed:realpart(ed) + %i*imagpart(ed), ed:ev(ed,numer), edp : -2*z*cos(z^2+c)+18*z^5+36*c*z^3+12*z^2+18*c^2*z+6*c*z+4*c, edp:expand(edp), edp:radcan(edp), edp:realpart(edp) + %i*imagpart(edp), edp:ev(edp,numer), z : z - (ed/edp), j : j+1 ), realpart(ed) + %i*imagpart(ed) )$ /* [wxMaxima: input end ] */ /* [wxMaxima: comment start ] [wxMaxima: comment end ] */ /* [wxMaxima: input start ] */ compile(iterdiffn)$ /* [wxMaxima: input end ] */ /* [wxMaxima: comment start ] This returns the value of the differential equation for the intial value of c. If this is close to zero then the last value of z might be a complex root of the Differential Equation . [wxMaxima: comment end ] */ /* [wxMaxima: input start ] */ n:6$ c:0.01 + %i*0.1$ iterdiffn(c,n); /* [wxMaxima: input end ] */ /* [wxMaxima: comment start ] [wxMaxima: comment end ] */ /* Maxima can't load/batch files which end with a comment! */ "Created with wxMaxima"$ Also the file iterdiff2c.mac /* . iterdiff2c.mac (c) Copyright 2017 , sciwise@ihug.co.nz Methodology for expressing differential equations in terms of an iterated function. In particular a first order differential equation and the derivative thereof . Note that f(z,c) might be any iterated function . . */ /* ' Differential equation we're examining . diff(f(z),z) + 3*f(z)*z = sin(z) for an iterated function z = f(z) , g(z) =f(f(z)) diff(g(z),z) + 3*g(z)*f(z) = sin(f(z)) . */ kill(all)$ eq1:diff(g(z,c),z) + 3*g(z,c)*f(z,c) - sin(f(z,c))$ expand(eq1); expand(diff(eq1,z)) ; /* For the iterated Mandelbrot fractal function . */ f(z,c):=z*z+c$ g(z,c):=f(f(z,c),c)$ f(z,c); g(z,c); /* */ eq1:diff(g(z,c),z) + 3*g(z,c)*f(z,c) - sin(f(z,c))$ ed:expand(eq1)$ edp:expand(diff(eq1,z))$ /* */ ed; edp; /* This may well simplify the construction of fractint formulas for this type of equation . */ /* z=Pixel c=Pixel : ed = -sin(z^2+c)+3*z^6+9*c*z^4+4*z^3+9*c^2*z^2+3*c*z^2+4*c*z+3*c^3+3*c^2 edp = -2*z*cos(z^2+c)+18*z^5+36*c*z^3+12*z^2+18*c^2*z+6*c*z+4*c */ kill(all)$ ratprint : false$ /* Iterated function , Differential and Newton's method. */ iterdiffn(c,n):=block([z,ed,edp,mag,j,i], mag:0.001, j:0, z:c, ed:1, for i:0 while( i < n) and ( abs(ed) > mag ) do( z:expand(z), z:radcan(z), z:realpart(z) + %i*imagpart(z), z:ev(z,numer), ed : -sin(z^2+c)+3*z^6+9*c*z^4+4*z^3+9*c^2*z^2+3*c*z^2+4*c*z+3*c^3+3*c^2, ed:expand(ed), ed:radcan(ed), ed:realpart(ed) + %i*imagpart(ed), ed:ev(ed,numer), edp : -2*z*cos(z^2+c)+18*z^5+36*c*z^3+12*z^2+18*c^2*z+6*c*z+4*c, edp:expand(edp), edp:radcan(edp), edp:realpart(edp) + %i*imagpart(edp), edp:ev(edp,numer), z : z - (ed/edp), j : j+1 ), realpart(ed) + %i*imagpart(ed) )$ /* */ compile(iterdiffn)$ /* This returns the value of the differential equation for the intial value of c. If this is close to zero then the last value of z might be a complex root of the Differential Equation . */ n:6$ c:0.01 + %i*0.1$ iterdiffn(c,n); /* */