User:Sam Derbyshire/MuPAD

Back to my user page, my talk page.

=MuPad code= Here is some code I used to generate some of my diagrams.

Be careful though, the code isn't exactly of high quality, it may not work as you'd expect it to.

General 2D function
f := plot::Function2d(airyAi(x), x=-6..4, LineWidth=0.8, LineColor=RGB::Red): g := plot::Function2d(airyBi(x), x=-6..4, LineWidth=0.8, LineColor=RGB::CornflowerBlue): plot(f,g, GridVisible = TRUE, SubgridVisible = TRUE, GridLineWidth = 0.2, SubgridLineWidth = 0.1, GridInFront = FALSE, GridLineColor = RGB::SlateGreyDark.[0.2], SubgridLineColor = RGB::SlateGreyDark.[0.1], TicksNumber = Low, TicksBetween = 4, TicksLabelFont=["Cambria",10,RGB::SlateGreyDark], AxesTitleFont=["Calibri",12,RGB::SlateGreyDark], ViewingBox = [-6..4, -2..8], Height=100, Width=100, Scaling=Constrained, AxesTitles = ["x", "y"] )

General 3D function
f := plot::Function3d(8*sin(x-cos(y))+(x^2+x*y),                    Mesh=[30,30], Submesh=[1,1]): plot(f,plot::Transform3d([0, 0, 0], [1, 0, 0, 0, 1, 0, 0, 0, 1], plot::modify(f, ZContours = [Automatic,14],                            LineWidth = 0.3,                             LineColorType = Dichromatic,                             LineColor = RGB::White.[0.2],                             LineColor2 = RGB::White.[0.2],                             XLinesVisible = FALSE,                             YLinesVisible = FALSE,                             Filled = FALSE)),       plot::Transform3d([0, 0, -15], [1, 0, 0, 0, 1, 0, 0, 0, 0], plot::modify(f, ZContours = [Automatic,14],                            LineWidth = 0.5,                             LineColorType = Dichromatic,                             LineColor = RGB::Red.[0.9],                             LineColor2 = RGB::CadmiumYellow.[0.9],                             XLinesVisible = FALSE,                             YLinesVisible = FALSE,                             Filled = FALSE)), GridVisible = TRUE, SubgridVisible = TRUE, GridLineWidth = 0.2, SubgridLineWidth = 0.1, FillColor=RGB::Red, FillColor2=RGB::CadmiumYellow, GridInFront = FALSE, GridLineColor = RGB::SlateGreyDark.[0.2], SubgridLineColor = RGB::SlateGreyDark.[0.1], TicksBetween = 1, TicksLabelFont=["Cambria",10,RGB::SlateGreyDark], AxesTitleFont=["Calibri",12,RGB::SlateGreyDark], Height=100, Width=100, ViewingBox=[-5..5,-5..5,-15..50], //Scaling = Constrained, AxesTitles = ["x", "y", "z"])

"[[Category:]]"==Complex conformal plots== m := 21: M := 2: f := z -> 1/z: hsv1 := z -> RGB::fromHSV([(arg(z))*180/PI, 0.5+abs(z)/(2*M*sqrt(2)), 0.25 + abs(z)/(4/3*M*sqrt(2)), 0.6]): hsv2 := z -> RGB::fromHSV([(arg(z))*180/PI, 0.5+abs(z)/(2*M*sqrt(2)), 0.25 + abs(z)/(4/3*M*sqrt(2)), 1]): p := arctan(Im(z)/Re(z)): q := abs(z)/(M*sqrt(2)): plot(plot::Conformal(z, z=-M*(1+I)..M*(1+I), Mesh=[m,m],LineColorType=Functional, LineColorFunction=(hsv1)), plot::Conformal(f(z), z=-M*(1+I)..M*(1+I), Mesh=[m,m], LineColorType=Functional, AdaptiveMesh=2, Submesh=[3,3], LineColorFunction=(hsv2)), Axes=Boxed, GridVisible = TRUE, SubgridVisible = TRUE,GridLineWidth = 0.1, SubgridLineWidth = 0.1, GridInFront = FALSE, GridLineColor = RGB::SlateGreyDark.[0.15], SubgridLineColor = RGB::SlateGreyDark.[0.1], ViewingBox = [-5..5, -5..5])

Contour plots
f := sin(PI*x)*sin(PI*y):

conts := 21: xlimit := 1: ylimit := 1: submeshlevel := 10:

a:=1: zmin := -1: zmax := 1: range:=zmax-zmin: plotmin := -1: plotmax := 1: projectionlevel:=plotmin: colour := x -> RGB::fromHSV([3*180/PI*x,1,1]): colfunc := (x,y,z) -> colour(floor((conts-1)/2*z)/((conts-1)/2)):

funcplot := plot::Function3d(f(x,y),                            x = -xlimit..xlimit,                             y = -ylimit..ylimit,                             Mesh = [21, 21],                             Submesh = [submeshlevel,submeshlevel],                             LineColor = RGB::Black.[0.4],                             LineWidth = 0.15,                             FillColorFunction = colfunc,                             AxesTitleFont=["Calibri",12,RGB::SlateGreyDark],                             ViewingBoxZRange = -1..1                            ):

contours := plot::modify(funcplot,                        ZContours = [Automatic, conts],                         LineWidth = 0.4,                         LineColor = RGB::Gray20.[0.9],                         XLinesVisible = FALSE,                         YLinesVisible = FALSE,                         Filled = FALSE                        ): ploteverything := plot::Canvas(funcplot, contours,                              Width = 200,                                Height = 200,                               TicksLabelFont=["Cambria",10,RGB::SlateGreyDark],                               plot::AmbientLight(1), Lighting=None, OrthogonalProjection = TRUE, plot::Camera([0, 0, 4*plotmax], [0, 0, (plotmin+plotmax)/2],PI/6)                              ): plot(ploteverything); Thanks to Inductiveload for helping out here.

Riemann Sphere
Simple projection of a curve (and coloured grid, cartesian/polar) back onto the Riemann sphere, can be done other way round. r := 1: V := 5: A := (X,Y) -> [X/(1+0.25*X^2+0.25*Y^2),Y/(1+0.25*X^2+0.25*Y^2),(-1+0.25*X^2+0.25*Y^2)/(1+0.25*X^2+0.25*Y^2)+1]: Am := (x,y,z) -> [2*x/(2-z), 2*y/(2-z)]: hsv1 := (t,x,y,z) -> RGB::fromHSV([(arg(x+I*y))*180/PI, 0.5+abs(x+I*y)/(2*V*sqrt(2)), 0.25 + abs(x+I*y)/(4/3*V*sqrt(2)), 0.9]): hsv2 := (t,x,y,z) -> RGB::fromHSV([(arg(Am(x,y,z)[1]+I*Am(x,y,z)[2]))*180/PI,    0.5+abs(Am(x,y,z)[1]+I*Am(x,y,z)[2])/(2*V*sqrt(2)), 0.25 + abs(Am(x,y,z)[1]+I*Am(x,y,z)[2])/(4/3*V*sqrt(2)), 0.9]): fgx := (i,j) -> (i,t,0)://fgx := (i,j) -> (i*cos(t),i*sin(t),0):// fgy := (i,j) -> (t,j,0)://fgy := (i,j) -> (cos(j*2*PI/V)*t,sin(j*2*PI/V)*t,0):// P := plot::Surface([u, v,0], u = -V .. V, v = -V .. V, UMesh=2*V+1, VMesh=2*V+1, Submesh=[1,1], FillColor=[1,1,1,1], FillColor2=[1,1,1,1], ULinesVisible=FALSE, VLinesVisible=FALSE): S := plot::Surface([cos(u)*sin(v), sin(u)*sin(v),cos(v)+1], u = 0 .. 2*PI, v = 0 .. PI,    Submesh=[1,1], FillColor=[1,1,1,1], FillColor2=[1,1,1,1], ULinesVisible=FALSE, VLinesVisible=FALSE): cx := t -> 1*(3*cos(t) - cos(3*t)): cy := t -> 1*(3*sin(t) - sin(3*t)): C := plot::Curve3d([cx(t),cy(t),0],t=0..2*PI, LineColor=RGB::Black.[0.9],LineWidth=0.8, Mesh=200): D := plot::Curve3d(A(cx(t),cy(t)),t=0..2*PI,LineColor=RGB::Black.[0.9],LineWidth=0.7, Mesh=100): Gx := i -> plot::Curve3d([fgx(i,j)],t=-V..V,LineWidth=0.4, LineColor=RGB::Black.[0.9], LineColorType=Functional, LineColorFunction=(hsv1)): Gy := j -> plot::Curve3d([fgy(i,j)],t=-V..V,LineWidth=0.4, LineColor=RGB::Black.[0.9], LineColorType=Functional, LineColorFunction=(hsv1)): Gsx := i -> plot::Curve3d(A(fgx(i,j)),t=-V..V,LineWidth=0.4, LineColor=RGB::Black.[0.9], LineColorType=Functional, LineColorFunction=(hsv2)): Gsy := j -> plot::Curve3d(A(fgy(i,j)),t=-V..V,LineWidth=0.4, LineColor=RGB::Black.[0.9], LineColorType=Functional, LineColorFunction=(hsv2)): Grid := plot::Group3d(Gx(k-V-1)$ k = 1..2*V+1 ,Gy(k-V-1) $ k = 1..2*V+1): SGrid := plot::Group3d(Gsx(k-V-1)$ k = 1..2*V+1 ,Gsy(k-V-1) $ k = 1..2*V+1): plot(P,S,C,D, Grid, SGrid,    ViewingBox=[-V..V,-V..V,-3..3], Height = 250, Width=250, Scaling=Constrained, AxesVisible=FALSE) Projection of the grid and the curve onto the Riemann Sphere, Euclidean motion of it ("rotation" matrix + translation) and projection back down. K := [-2,0,3]: R := [K[1],K[2],K[3]+1]: tra := [K[1],K[2],K[3]]: V := 6: T := matrix(0.707,0,0.707],[0,1,0],[-0.707,0,0.707): p1 := plot::Point3d(K, PointSize = 2, PointColor=RGB::Black): p2 := plot::Point3d(R, PointSize = 2, PointColor=RGB::Black): Ab := (X,Y) -> [2*X/(1+X^2+Y^2),2*Y/(1+X^2+Y^2),(-1+X^2+Y^2)/(1+X^2+Y^2)]: Af := (x,y,z) -> [x/(1.01-z), y/(1.01-z)]: Ab1 := (X,Y) -> Ab(X,Y): Af1 := (x,y,z) -> Af(x,y,z): Ab2 := (X,Y) -> zip(Ab((X-K[1])/R[3],(Y-K[2])/R[3]),[K[1],K[2],K[3]], _plus): Af2 := (x,y,z) -> zip(Af(R[3]*(x-K[1]),R[3]*(y-K[2]),z-K[3]),[K[1],K[2]], _plus): hsv := (t,x,y,z) -> RGB::fromHSV([(arg(x+I*y))*180/PI, 0.5+abs(x+I*y)/(2*V*sqrt(2)), 0.25 + abs(x+I*y)/(4/3*V*sqrt(2)), 0.3]): hsv1 := (t,x,y,z) -> RGB::fromHSV([(arg(x+I*y))*180/PI, 0.5+abs(x+I*y)/(2*V*sqrt(2)), 0.25 + abs(x+I*y)/(4/3*V*sqrt(2)), 0.9]): hsv2 := (t,x,y,z) -> hsv1(t,Af1(x,y,z)[1],Af1(x,y,z)[2],z): hsv3 := (t,x,y,z) -> hsv2(t,(T^(-1)*(matrix([Ab2(x,y)[1],Ab2(x,y)[2],Ab2(x,y)[3]])-matrix(tra)))[1],    (T^(-1)*(matrix([Ab2(x,y)[1],Ab2(x,y)[2],Ab2(x,y)[3]])-matrix(tra)))[2],     (T^(-1)*(matrix([Ab2(x,y)[1],Ab2(x,y)[2],Ab2(x,y)[3]])-matrix(tra)))[3]): fgx := (i,j) -> (i*cos(t),i*sin(t),0)://fgx := (i,j) -> (i,t,0):// fgy := (i,j) -> (cos(j*2*PI/V)*t,sin(j*2*PI/V)*t,0)://fgy := (i,j) -> (t,j,0):// P := plot::Surface([u, v,0], u = -V .. V, v = -V .. V, UMesh=2*V+1, VMesh=2*V+1, Submesh=[1,1], FillColor=[1,1,1,1], FillColor2=[1,1,1,1], ULinesVisible=FALSE, VLinesVisible=FALSE): S1 := plot::Surface([cos(u)*sin(v), sin(u)*sin(v),cos(v)], u = 0 .. 2*PI, v = 0 .. PI,    Submesh=[1,1], FillColor=[1,1,1,1], FillColor2=[1,1,1,1], ULinesVisible=FALSE, VLinesVisible=FALSE): S2 := plot::Surface([cos(u)*sin(v)+K[1], sin(u)*sin(v)+K[2],cos(v)+K[3]], u = 0 .. 2*PI, v = 0 .. PI,    Submesh=[1,1], FillColor=[1,1,1,1], FillColor2=[1,1,1,1], ULinesVisible=FALSE, VLinesVisible=FALSE): cx := t -> 2.25*cos(t)-0.7*cos(13*t): cy := t -> 2.25*sin(t)-0.7*sin(13*t): C1 := plot::Curve3d([cx(t),cy(t),0],t=0..2*PI, LineColor=RGB::Black.[0.3],LineWidth=0.8, Mesh=400): CS := plot::Curve3d(Ab1(cx(t),cy(t)),t=0..2*PI,LineColor=RGB::Black.[0.9],LineWidth=0.7, Mesh=500): CS2 := plot::Transform3d(tra,T,CS): C2 := plot::Curve3d(Af2((T*matrix([Ab1(cx(t),cy(t))[1],Ab1(cx(t),cy(t))[2],Ab1(cx(t),cy(t))[3]])+matrix(tra))[1], (T*matrix([Ab1(cx(t),cy(t))[1],Ab1(cx(t),cy(t))[2],Ab1(cx(t),cy(t))[3]])+matrix(tra))[2], (T*matrix([Ab1(cx(t),cy(t))[1],Ab1(cx(t),cy(t))[2],Ab1(cx(t),cy(t))[3]])+matrix(tra))[3]).[0],    t=0..2*PI, LineColor=RGB::Black.[0.9],LineWidth=0.8, Mesh=600): Gx := i -> plot::Curve3d([fgx(i,j)],t=-V..V,LineWidth=0.4, LineColor=RGB::Black.[0.9], LineColorType=Functional, LineColorFunction=(hsv)): Gy := j -> plot::Curve3d([fgy(i,j)],t=-V..V,LineWidth=0.4, LineColor=RGB::Black.[0.9], LineColorType=Functional, LineColorFunction=(hsv)): Gsx := i -> plot::Curve3d(Ab1(fgx(i,j)),t=-V..V,LineWidth=0.4, LineColor=RGB::Black.[0.9], LineColorType=Functional, LineColorFunction=(hsv2)): Gsy := j -> plot::Curve3d(Ab1(fgy(i,j)),t=-V..V,LineWidth=0.4, LineColor=RGB::Black.[0.9], LineColorType=Functional, LineColorFunction=(hsv2)): G2x := i -> plot::Curve3d(Af2((T*matrix(Ab1(fgx(i,j)))+matrix(tra))[1],(T*matrix(Ab1(fgx(i,j)))+matrix(tra))[2],(T*matrix(Ab1(fgx(i,j)))+matrix(tra))[3]).[0],    t=-V..V,LineWidth=0.4, LineColor=RGB::Black.[0.9], LineColor=RGB::Black.[0.9], LineColorType=Functional, LineColorFunction=(hsv3)): G2y := j -> plot::Curve3d(Af2((T*matrix(Ab1(fgy(i,j)))+matrix(tra))[1],(T*matrix(Ab1(fgy(i,j)))+matrix(tra))[2],(T*matrix(Ab1(fgy(i,j)))+matrix(tra))[3]).[0],    t=-V..V,LineWidth=0.4, LineColor=RGB::Black.[0.9], LineColor=RGB::Black.[0.9], LineColorType=Functional, LineColorFunction=(hsv3)): Grid := plot::Group3d(Gx(k-V-1)$ k = 1..2*V+1 ,Gy(k-V-1) $ k = 1..2*V+1): SGrid := plot::Group3d(Gsx(k-V-1)$ k = 1..2*V+1 ,Gsy(k-V-1) $ k = 1..2*V+1): SGrid2 := plot::Group3d(plot::Transform3d(tra,T,Gsx(k-V-1))$ k = 1..2*V+1 ,plot::Transform3d(tra,T,Gsy(k-V-1)) $ k = 1..2*V+1): Grid2 := plot::Group3d(G2x(k-V-1)$ k = 1..2*V+1 ,G2y(k-V-1) $ k = 1..2*V+1): plot(P,S2,C1,CS2,C2, Grid, Grid2, SGrid2, p1,p2, S1,SGrid,CS,    ViewingBox=[-V-3..V+3,-V-3..V+3,-4..4], Height = 250, Width=250, Scaling=Constrained, AxesVisible=FALSE)

Epitrochoids and Hypotrochoids
r := 2.1: R := 3: s := -1: // 1 for epi, -1 for hypo displacement := 0: angle := 0: period := 14*PI: range := R + 2*r*heaviside(s)+displacement*heaviside(displacement): centerx := a -> (R+s*r)*cos(+a): centery := a -> (R+s*r)*sin(+a): center2x := a -> cos(angle)*centerx(a) - sin(angle)*centery(a): center2y := a -> sin(angle)*centerx(a) + cos(angle)*centery(a): pointx := t -> center2x(t)+(r+displacement)*cos(s*(R/r+s*1)*t): pointy := t -> center2y(t)+(r+displacement)*sin(s*(R/r+s*1)*t): point2x := t -> cos(-angle)*pointx(t) - sin(-angle)*pointy(t): point2y := t -> sin(-angle)*pointx(t) + cos(-angle)*pointy(t):

plot( plot::Circle2d(R, LineColor=RGB::Blue, LineWidth=0.6), plot::Circle2d(r, [centerx(a),centery(a)], a=0..period, LineColor=RGB::Black, LineWidth=0.6), plot::Curve2d([point2x(t),point2y(t)], t=0..a, a=0..period, LineWidth=0.6, LineColor=RGB::Red, Mesh=500), plot::Line2d([centerx(a),centery(a)],[point2x(a),point2y(a)], a=0..period, LineColor=RGB::Black), plot::Point2d([centerx(a),centery(a)], a=0..period, PointSize=2, PointColor=RGB::Black), plot::Point2d([point2x(a),point2y(a)], a=0..period, PointSize=2, PointColor=RGB::Red),

GridVisible = TRUE, SubgridVisible = TRUE, GridLineWidth = 0.2, SubgridLineWidth = 0.1, GridInFront = FALSE, GridLineColor = RGB::SlateGreyDark.[0.2], SubgridLineColor = RGB::SlateGreyDark.[0.1], TicksNumber = Low, TicksBetween = 1, TicksLabelFont=["Cambria",10,RGB::SlateGreyDark], AxesTitleFont=["Calibri",12,RGB::SlateGreyDark], ViewingBox = [-range..range, -range..range], Frames=100, TimeEnd=5, Height=200, Width=200, Scaling=Constrained, AxesTitles = ["x", "y"] )

Pedal curve
fx := t -> 2*cos(t): fy := t -> sin(t): f := t -> [fx(t),fy(t)]: range :=2.5: px := 0: py := 0: h := 0.1: slope := a -> piecewise([diff(fx(a),a)=0,10^10],[diff(fx(a),a) <> 0,diff(fy(a),a)/(diff(fx(a),a))]): angle := a -> arctan(slope(a)): ix := a ->((slope(a)*(-fx(a))+fy(a))-(-(1/slope(a))*(-px)+py))/(-slope(a)-1/slope(a)): iy := a -> slope(a)*(ix(a)-fx(a))+fy(a): period := 2*PI: plot( plot::Curve2d(f(t),t=0..2*PI, LineColor=RGB::Black, LineWidth=0.6, Mesh=200), plot::Function2d(slope(a)*(x-fx(a))+fy(a),x=-range..range, a=0..period, LineColor=RGB::Blue, LineWidth=0.5), plot::Function2d(-(1/slope(a))*(x-px)+py,x=-range..range, a=0..period, LineColor=RGB::SapGreen, LineWidth=0.5), plot::Curve2d([ix(t),iy(t)],t=0..period, LineWidth=0.6, LineColorFunction = ((u, x, y, a) -> [1, 0, 0,1]), Mesh=400, Submesh=3, AdaptiveMesh=2), plot::Point2d([px,py], PointSize=2, PointColor=RGB::Black), plot::Point2d([fx(a),fy(a)], a=0..period, PointSize=2, PointColor=RGB::Black), plot::Point2d([ix(a),iy(a)], a=0..period, PointSize=2, PointColor=RGB::Red), plot::Polygon2d(ix(a)+h*cos(angle(a)), iy(a)+h*cos(angle(a))*slope(a)],                [ix(a)+h*cos(angle(a))+h*sin(angle(a)), iy(a)+h*cos(angle(a))*slope(a)-h*sin(angle(a))/slope(a)],                 [ix(a)+h*sin(angle(a)), iy(a)-h*sin(angle(a))/slope(a), a=0..period, LineColor=RGB::Black, LineWidth=0.5), GridVisible = TRUE, SubgridVisible = TRUE, GridLineWidth = 0.2, SubgridLineWidth = 0.1, GridInFront = FALSE, GridLineColor = RGB::SlateGreyDark.[0.2], SubgridLineColor = RGB::SlateGreyDark.[0.1], TicksNumber = Low, TicksBetween = 1, TicksLabelFont=["Cambria",10,RGB::SlateGreyDark], AxesTitleFont=["Calibri",12,RGB::SlateGreyDark], ViewingBox = [-range..range, -range..range], Frames=400, TimeEnd=10, Height=200, Width=200, Scaling=Constrained, AxesTitles = ["x", "y"] )

Orthotomic curve
fx := t -> 2*cos(t): fy := t -> sin(t): f := t -> [fx(t),fy(t)]: range :=4.5: px := 0: py := 1: h := 0.1: slope := a -> piecewise([diff(fx(a),a)=0,10^10],[diff(fx(a),a) <> 0,diff(fy(a),a)/(diff(fx(a),a))]): angle := a -> arctan(slope(a)): ix := a ->((slope(a)*(-fx(a))+fy(a))-(-(1/slope(a))*(-px)+py))/(-slope(a)-1/slope(a)): iy := a -> slope(a)*(ix(a)-fx(a))+fy(a): period := 2*PI: plot( plot::Curve2d(f(t),t=0..2*PI, LineColor=RGB::Black, LineWidth=0.6, Mesh=200), plot::Function2d(slope(a)*(x-fx(a))+fy(a),x=-range..range, a=0..period, LineColor=RGB::Blue, LineWidth=0.5), plot::Function2d(-(1/slope(a))*(x-px)+py,x=-range..range, a=0..period, LineColor=RGB::SapGreen, LineWidth=0.5), plot::Curve2d([px+2*(ix(t)-px),py+2*(iy(t)-py)],t=0..period, LineWidth=0.6, LineColorFunction = ((u, x, y, a) -> [1, 0, 0,1]), Mesh=400, Submesh=3, AdaptiveMesh=2), plot::Point2d([px,py], PointSize=2, PointColor=RGB::Black), plot::Point2d([fx(a),fy(a)], a=0..period, PointSize=2, PointColor=RGB::Black), plot::Point2d([ix(a),iy(a)], a=0..period, PointSize=2, PointColor=RGB::Blue), plot::Point2d([px+2*(ix(a)-px),py+2*(iy(a)-py)], a=0..period, PointSize=2, PointColor=RGB::Red), plot::Polygon2d(ix(a)+h*cos(angle(a)), iy(a)+h*cos(angle(a))*slope(a)],                [ix(a)+h*cos(angle(a))+h*sin(angle(a)), iy(a)+h*cos(angle(a))*slope(a)-h*sin(angle(a))/slope(a)],                 [ix(a)+h*sin(angle(a)), iy(a)-h*sin(angle(a))/slope(a), a=0..period, LineColor=RGB::Black, LineWidth=0.5), GridVisible = TRUE, SubgridVisible = TRUE, GridLineWidth = 0.2, SubgridLineWidth = 0.1, GridInFront = FALSE, GridLineColor = RGB::SlateGreyDark.[0.2], SubgridLineColor = RGB::SlateGreyDark.[0.1], TicksNumber = Low, TicksBetween = 1, TicksLabelFont=["Cambria",10,RGB::SlateGreyDark], AxesTitleFont=["Calibri",12,RGB::SlateGreyDark], ViewingBox = [-range..range, -range..range], Frames=400, TimeEnd=10, Height=200, Width=200, Scaling=Constrained, AxesTitles = ["x", "y"] )

Contrapedal curve
fx := t -> 2*cos(t): fy := t -> sin(t): f := t -> [fx(t),fy(t)]: range :=2.5: px := 0: py := 0: h := 0.1: slope := a -> piecewise([diff(fy(a),a)=0,-10^10],[diff(fy(a),a) <> 0,-diff(fx(a),a)/(diff(fy(a),a))]): angle := a -> arctan(slope(a)): ix := a ->((slope(a)*(-fx(a))+fy(a))-(-(1/slope(a))*(-px)+py))/(-slope(a)-1/slope(a)): iy := a -> slope(a)*(ix(a)-fx(a))+fy(a): period := 2*PI: plot( plot::Curve2d(f(t),t=0..2*PI, LineColor=RGB::Black, LineWidth=0.6, Mesh=200), plot::Function2d(slope(a)*(x-fx(a))+fy(a),x=-range..range, a=0..period, LineColor=RGB::Blue, LineWidth=0.5), plot::Function2d(-(1/slope(a))*(x-px)+py,x=-range..range, a=0..period, LineColor=RGB::SapGreen, LineWidth=0.5), plot::Curve2d([ix(t),iy(t)],t=0..period, LineWidth=0.6, LineColorFunction = ((u, x, y, a) -> [1, 0, 0,1]), Mesh=400, Submesh=3, AdaptiveMesh=2), plot::Point2d([fx(a),fy(a)], a=0..period, PointSize=2, PointColor=RGB::Black), plot::Point2d([px,py], PointSize=2, PointColor=RGB::Black), plot::Point2d([ix(a),iy(a)], a=0..period, PointSize=2, PointColor=RGB::Red), plot::Polygon2d(ix(a)+h*cos(angle(a)), iy(a)+h*cos(angle(a))*slope(a)],                [ix(a)+h*cos(angle(a))+h*sin(angle(a)), iy(a)+h*cos(angle(a))*slope(a)-h*sin(angle(a))/slope(a)],                 [ix(a)+h*sin(angle(a)), iy(a)-h*sin(angle(a))/slope(a), a=0..period, LineColor=RGB::Black, LineWidth=0.5), GridVisible = TRUE, SubgridVisible = TRUE, GridLineWidth = 0.2, SubgridLineWidth = 0.1, GridInFront = FALSE, GridLineColor = RGB::SlateGreyDark.[0.2], SubgridLineColor = RGB::SlateGreyDark.[0.1], TicksNumber = Low, TicksBetween = 1, TicksLabelFont=["Cambria",10,RGB::SlateGreyDark], AxesTitleFont=["Calibri",12,RGB::SlateGreyDark], ViewingBox = [-range..range, -range..range], Frames=400, TimeEnd=10, Height=200, Width=200, Scaling=Constrained, AxesTitles = ["x", "y"] )

Evolute
fx := t -> 2*cos(t): fy := t -> sin(t): f := t -> [fx(t),fy(t)]: range :=5: slope := a -> piecewise([diff(fx(a),a)=0,10^10],[diff(fx(a),a) <> 0,diff(fy(a),a)/(diff(fx(a),a))]): angle := a -> arctan(slope(a)): period := 2*PI: lines_number := 50: timend := 5:

lines := [FAIL $ lines_number+1]: for i from 0 to lines_number do   theta := i*2*PI/(lines_number): lines[i+1] := plot::Function2d((-1/slope(t)*(x-fx(t))+fy(t))|(t=theta),                   VisibleAfter = timend*i/(lines_number+1),                   Color = RGB::Blue                     ): end_for:

cx :=t -> fx(t)+diff(fy(t),t)*((diff(fx(t),t))^2+(diff(fy(t),t))^2)/(diff(fx(t),t,t)*diff(fy(t),t)-diff(fy(t),t,t)*diff(fx(t),t)): cy :=t -> fy(t)+diff(fx(t),t)*((diff(fx(t),t))^2+(diff(fy(t),t))^2)/(diff(fy(t),t,t)*diff(fx(t),t)-diff(fx(t),t,t)*diff(fy(t),t)): plot( plot::Curve2d(f(t),t=0..2*PI, LineColor=RGB::Black, LineWidth=0.6, Mesh=200), plot::Group2d(op(lines)), plot::Curve2d([cx(t),cy(t)],t=0..a,a=0..2*PI, LineColor=RGB::Blue, LineWidth=0.6, Mesh=200), plot::Point2d([cx(a),cy(a)],a=0..2*PI, PointSize=2,Color=RGB::Black), plot::Point2d([fx(a),fy(a)], a=0..period, PointSize=2, PointColor=RGB::Black), GridVisible = TRUE, SubgridVisible = TRUE, GridLineWidth = 0.2, SubgridLineWidth = 0.1, GridInFront = FALSE, GridLineColor = RGB::SlateGreyDark.[0.2], SubgridLineColor = RGB::SlateGreyDark.[0.1], TicksNumber = Low, TicksBetween = 1, TicksLabelFont=["Cambria",10,RGB::SlateGreyDark], AxesTitleFont=["Calibri",12,RGB::SlateGreyDark], ViewingBox = [-range..range, -range..range], Frames=400, TimeEnd=timend, Height=200, Width=200, Scaling=Constrained, AxesTitles = ["x", "y"] )

With parallel curves
fx := t -> 2*cos(t): fy := t -> sin(t): f := t -> [fx(t),fy(t)]: g := (t,a) -> [fx(t) + a*diff(fy(t),t)/(sqrt((diff(fx(t),t))^2+(diff(fy(t),t)^2))),fy(t) - a*diff(fx(t),t)/(sqrt((diff(fx(t),t))^2+(diff(fy(t),t)^2)))]: range :=5: slope := a -> piecewise([diff(fx(a),a)=0,10^10],[diff(fx(a),a) <> 0,diff(fy(a),a)/(diff(fx(a),a))]): angle := a -> arctan(slope(a)): period := 2*PI: lines_number := 21: timend := 5:

lines := [FAIL $ lines_number+1]: for i from 0 to lines_number do   theta := i*2*PI/(lines_number): lines[i+1] := plot::Function2d((-1/slope(t)*(x-fx(t))+fy(t))|(t=theta),                   VisibleAfter = timend*i/(lines_number+1),                   Color = RGB::Blue,                   TimeEnd = 5                     ): end_for:

cx :=t -> fx(t)+diff(fy(t),t)*((diff(fx(t),t))^2+(diff(fy(t),t))^2)/(diff(fx(t),t,t)*diff(fy(t),t)-diff(fy(t),t,t)*diff(fx(t),t)): cy :=t -> fy(t)+diff(fx(t),t)*((diff(fx(t),t))^2+(diff(fy(t),t))^2)/(diff(fy(t),t,t)*diff(fx(t),t)-diff(fx(t),t,t)*diff(fy(t),t)): plot( plot::Curve2d(f(t),t=0..2*PI, LineColor=RGB::Red, LineWidth=0.6, Mesh=200), plot::Group2d(op(lines)), plot::Curve2d(g(t,a),t=0..2*PI,a=0..-5, LineColor=RGB::Black, LineWidth=0.6, Mesh=200, TimeBegin=5, TimeEnd=10, VisibleBeforeBegin = FALSE), //plot::Curve2d([cx(t),cy(t)],t=0..a,a=0..2*PI, LineColor=RGB::Blue, LineWidth=0.6, Mesh=200), //plot::Point2d([cx(a),cy(a)],a=0..2*PI, PointSize=2,Color=RGB::Black), plot::Point2d([fx(a),fy(a)], a=0..period, PointSize=2, PointColor=RGB::Black), GridVisible = TRUE, SubgridVisible = TRUE, GridLineWidth = 0.2, SubgridLineWidth = 0.1, GridInFront = FALSE, GridLineColor = RGB::SlateGreyDark.[0.2], SubgridLineColor = RGB::SlateGreyDark.[0.1], TicksNumber = Low, TicksBetween = 1, TicksLabelFont=["Cambria",10,RGB::SlateGreyDark], AxesTitleFont=["Calibri",12,RGB::SlateGreyDark], ViewingBox = [-range..range, -range..range], Frames=400, TimeEnd=timend, Height=200, Width=200, Scaling=Constrained, AxesTitles = ["x", "y"] )

Arclength parametrization case
Disclaimer : Works only with arclenth parametrization ! fx := t -> arcsinh(t): fy := t -> sqrt(t^2+1): f := t -> [fx(t),fy(t)]: range :=4: slope := a -> piecewise([diff(fx(a),a)=0,10^10],[diff(fy(a),a)=0,10^-10],[diff(fx(a),a) <> 0 and diff(fy(a),a) <> 0,diff(fy(a),a)/(diff(fx(a),a))]): angle := a -> arctan(slope(a)): period := 2*PI: lines_number := 50: timend := 5: cx := t -> fx(t)-t*diff(fx(t),t): cy := t -> fy(t)-t*diff(fy(t),t): plot( plot::Curve2d(f(t),t=-6..6, LineColor=RGB::Blue, LineWidth=0.6, Mesh=200), plot::Point2d([cx(a),cy(a)], a=-5..5, PointSize=2, PointColor=RGB::Black, VisibleFromTo=0..timend), plot::Line2d([fx(a),fy(a)],[cx(a),cy(a)], a=-5..5, LineWidth=0.4, LineColor=RGB::Black, VisibleFromTo=0..timend), plot::Curve2d([cx(t),cy(t)], t=-30..30, LineWidth=0.6, LineColor=[1,0.8,0.8], Mesh=200), plot::Curve2d([cx(t),cy(t)], t=-5..a,a=-5..5, LineWidth=0.6, LineColor=RGB::Red, Mesh=200), plot::Curve2d([cx(t),cy(t)], t=-30..30, LineWidth=0.6, LineColor=RGB::Red, VisibleAfter=timend, Mesh=200), GridVisible = TRUE, SubgridVisible = TRUE, GridLineWidth = 0.2, SubgridLineWidth = 0.1, GridInFront = FALSE, GridLineColor = RGB::SlateGreyDark.[0.2], SubgridLineColor = RGB::SlateGreyDark.[0.1], TicksNumber = Low, TicksBetween = 1, TicksLabelFont=["Cambria",10,RGB::SlateGreyDark], AxesTitleFont=["Calibri",12,RGB::SlateGreyDark], ViewingBox = [-3..3, -1..5], Frames=400, TimeEnd=timend, Height=160, Width=120, Scaling=Constrained, AxesTitles = ["x", "y"] )

Variant, all parametrizations
To be used against an evolute : first equations are those of the evolute. fx := t -> (3+1)*(cos(t)-cos((3+1)*t)/(3+1)): fy := t -> (3+1)*(sin(t)-sin((3+1)*t)/(3+1)): f := t -> [fx(t),fy(t)]: range :=4: slope := a -> piecewise([diff(fx(a),a)=0,10^10],[diff(fy(a),a)=0,10^-10],[diff(fx(a),a) <> 0 and diff(fy(a),a) <> 0,diff(fy(a),a)/(diff(fx(a),a))]): angle := a -> arctan(slope(a)): period := 2*PI: lines_number := 100: timend := 10: cx :=t -> fx(t)+diff(fy(t),t)*((diff(fx(t),t))^2+(diff(fy(t),t))^2)/(diff(fx(t),t,t)*diff(fy(t),t)-diff(fy(t),t,t)*diff(fx(t),t)): cy :=t -> fy(t)+diff(fx(t),t)*((diff(fx(t),t))^2+(diff(fy(t),t))^2)/(diff(fy(t),t,t)*diff(fx(t),t)-diff(fx(t),t,t)*diff(fy(t),t)): plot( plot::Curve2d(f(t),t=-PI..a,a=-PI..PI, LineColor=RGB::Red, LineWidth=0.6, Mesh=200), //plot::Curve2d([-0.6*fx(t),-0.6*fy(t)],t=-2*PI..2*PI, LineColor=RGB::Purple, LineWidth=0.6, Mesh=200), plot::Curve2d([cx(t),cy(t)],t=-2*PI..2*PI, LineColor=RGB::Blue, LineWidth=0.6, Mesh=200), plot::Point2d([fx(a),fy(a)], a=-PI..PI, PointSize=2, PointColor=RGB::Black), //plot::Point2d([cx(a),cy(a)], a=-PI..PI, PointSize=2, PointColor=RGB::Black), plot::Line2d([cx(a),cy(a)],[fx(a),fy(a)], a=-PI..PI, LineColor=RGB::Black, LineWidth=0.4), GridVisible = TRUE, SubgridVisible = TRUE, GridLineWidth = 0.2, SubgridLineWidth = 0.1, GridInFront = FALSE, GridLineColor = RGB::SlateGreyDark.[0.2], SubgridLineColor = RGB::SlateGreyDark.[0.1], TicksNumber = Low, TicksBetween = 1, TicksLabelFont=["Cambria",10,RGB::SlateGreyDark], AxesTitleFont=["Calibri",12,RGB::SlateGreyDark], ViewingBox = [-6..6, -6..6], Frames=400, TimeEnd=timend, Height=160, Width=120, Scaling=Constrained, AxesTitles = ["x", "y"] )

Two coupled nonlinear differential equations
Here Lotka-Volterra equations with decreasing returns. [r0,r1,r2,r3,r4,r5] := [10, 4.00000, 0.00000,-0.50000,-0.05000, 0.00000]: [h0,h1,h2,h3,h4,h5] := [10,-1.00000, 0.00000, 0.30000, 0.00000, 0.00000]: fr := (yr,yh) ->            r1*yr + r2*yh + r3*yr*yh + r4*yr^2 + r5*yh^2: fh := (yr,yh) ->            h1*yh + h2*yr + h3*yh*yr + h4*yh^2 + h5*yr^2: Y := [yr,yh]: Y0 := [r0,h0]: Z[0]:=Y0: f := (t,Y) -> [fr(Y[1],Y[2]),fh(Y[1],Y[2])]:

N := 800:

for i from 0 to N do  t[i] := i/20 end_for:

for i from 1 to N do   Z[i] := numeric::odesolve(f, t[i-1]..t[i], Z[i-1]) end_for:

plot(plot::Listplot([ [ t[i],Z[i][1] ] $ i = 0..N ], LineWidth=0.8, LineColor=RGB::SapGreen), plot::Listplot([ [ t[i],Z[i][2] ] $ i = 0..N ], LineWidth=0.8, LineColor=RGB::Red), PointsVisible=FALSE, AxesTitles = ["Time", "Population"], GridVisible = TRUE, SubgridVisible = TRUE, GridLineWidth = 0.2, SubgridLineWidth = 0.1, GridInFront = FALSE, GridLineColor = RGB::SlateGreyDark.[0.2], SubgridLineColor = RGB::SlateGreyDark.[0.1]):

plot(plot::Listplot([ [ Z[i][1],Z[i][2] ] $ i = 0..N ], LineWidth=0.6,LineColor=RGB::Blue), PointsVisible=FALSE, AxesTitles = ["Prey", "Predator"], GridVisible = TRUE, SubgridVisible = TRUE, GridLineWidth = 0.2, SubgridLineWidth = 0.1, GridInFront = FALSE, GridLineColor = RGB::SlateGreyDark.[0.2], SubgridLineColor = RGB::SlateGreyDark.[0.1]):