Second order non-linear equation solver in Commodore BASIC v3.5

diff --git a/solver.bas b/solver.bas
new file mode 100644 (file)
index 0000000..8029556
--- /dev/null
+++ b/solver.bas
@@ -0,0 +1,84 @@
+10 REM Y.H.KU NON-LINEAR EQUATION SOLVER
+20 REM ENTER EQUATION AT LINE 800
+30 COLOR 1,1:COLOR 0,2,4:COLOR 4,1:SCNCLR
+40 OPEN3,6:OPEN 5,3,3:OPEN2,6,2:OPEN1,6,1
+45 PRINT#1,"H"
+50 PRINT"***************************************"
+60 PRINT"*                                     *"
+70 PRINT"*NON-LINEAR 2ND. ORDER EQUATION SOLVER*"
+80 PRINT"*                                     *"
+90 PRINT"*   USING Y.H.KU'S ACC. PHASE PLANE   *"
+100 PRINT"*                                     *"
+110 PRINT"*  M.K.HOBDEN,PLUS-4 VERSION MAY 89   *"
+120 PRINT"*                                     *"
+130 PRINT"***************************************"
+135 PRINT
+140 PRINT"  ENTER D.E AT LINE 800."
+150 PRINT
+160 PRINT" ENTER INITIAL VELOCITY AND DISPLACEMENT"
+170 PRINT" AT THE PROGRAM PROMPTS"
+180 PRINT
+190 INPUT" INITIAL DISP.(X0)";X0
+200 PRINT
+210 INPUT" INITIAL VEL. (V0)";V0
+220 PRINT
+230 INPUT" CONT. PARAM. (MU)";MU
+240 PRINT
+245 INPUT" STEP SIZE    (<1)";D
+250 DIM A(300),E(300)
+260 PRINT:PRINT"      CALCULATING......"
+270 IF V0<0 THEN 320
+280 IF V0>0 THEN 300
+290 IF X0>0 THEN 320
+300 Z=D
+310 GOTO 330
+320 Z=-D
+330 X1=X0+Z
+340 GOSUB 800
+350 C=Z*X1
+360 S=(B*B)-4*C
+370 IF S<0 THEN 540
+380 K=Z/ABS(Z)
+390 V1=(-B+K*SQR((B*B)-4*C))/2
+400 IF G%=1 THEN 430
+410 IF INT(H/50)<>H/50 THEN 430
+430 I=I+1
+440 A(I)=X1
+450 E(I)=V1
+460 H=H+1
+470 IF H=50 THEN 590
+480 V0=V1
+490 X0=X1
+500 IF V0=0 THEN 520
+510 GOTO 330
+520 IF E(I-1)<0 THEN 300
+530 GOTO 320
+540 Z=Z/10
+550 IF ABS(Z)<1.0E-3 THEN 570
+560 GOTO 330
+570 V1=0
+580 GOTO 430
+590 IF J>0 THEN 660
+600 SCNCLR:PRINT:PRINT
+610 PRINT"   SOLUTION OF EQUATION:"
+620 PRINT
+630 PRINT"   ";Q$
+640 PRINT:PRINT"    FOR MU=";MU:PRINT
+650 PRINT:PRINT"    USING Y.H.KU'S METHOD"
+660 IF G%=1 THEN 720
+670 PRINT#1,"D",479,0
+680 PRINT#1,"M",240,-240:PRINT#1,"D",240,240
+690 PRINT#1,"M",0,240:PRINT#3,"PHASE PLANE SOLUTION":PRINT#1,"M",245,0:PRINT#3,"X'"
+700 PRINT#1,"M",470,-180:PRINT#3,"X"
+710 PRINT#1,"M",240,0:PRINT#1,"I"
+720 IF G%=0 THEN GOSUB 1000
+730 FOR I=1 TO H
+740 PRINT#1,"J",0+(A(I)*-75),(E(I)*-75)
+750 NEXT I
+760 H=0:I=0:J=J+1:GOTO 480
+800 B=-(V0+Z*(MU*(1-(X1*X1))))
+810 Q$="X#-MU(1-(X*X))X'+X=0"
+820 RETURN
+1000 PRINT#1,"R",0+(A(1)*-75),(E(1)*-75)
+1100 G%=1:RETURN
+1200 END