For those of you (if any) expecting a solution using just some tiny HP calculator, read on, I include a comprehensive one at the end of this post, but first some pertinent comments:

dbrunell posted:

*"The solution to the equation can be represented by: [an explicit symbolic solution found using Mathematica]"*

First of all, let me clearly state that I thank you very much for your interest in this informal 'challenge' I issued, and rest assured that none of my remarks are meant to be derogatory of your efforts to solve it. It's just that I am the kind of person that likes to abide by the rules of the game and expect others to do likewise, so I'm somewhat upset when someone disregards them, specially when it's totally unnecessary.

As you may guess, using Mathematica, MathLab, Derive, or whatever to solve one of my challenges is *utterly useless* as far as showing and learning new tricks with out beloved HP calculators is concerned, and as far as I'm concerned as well. As I said before, I also own a copy of Mathematica, and can easily achieve these results with no effort and to no purpose whatsoever.

Should you had found the explicit symbolic solution using some CAS or symbolic package in your HP48/49, that would be a thing. But finding it using Mathematica is breaking the rules big time. Now, once Mathematica has given you that weird-looking expression in terms of Bessel functions , you then go on:

*"Here is my 48S program for Bessel function of the first kind: "*

**THIS** is more like it !! That's what I want, some real code running in a real (or emulated) HP calculator, solving some difficult problem. What a pity that you chose to resort to that Mathematica of yours, instead of fully solving it with your powerful 48S from the scratch !

But you then end it all saying:

*"I will leave it to you to write a Bessel function for the HP-25. :) "*

What for ? WHO needs Bessel functions to solve this problem ? Mathematica needs them ? And WHO needs Mathematica for this problem ? Not anyone owning most any HP RPN or RPL calculators, for certain ...

In order to prove my point, here's a complete, reasoned solution to this challenge, using nothing more and nothing else than a vintage HP calculator. The solution is given for the HP-15C, for the one and only reason that it's the HP model I use routinely. A similar solution can be given for most other models, up to *and including* the HP-25 itself, but I don't have an HP-25 at hand right now. And yes, you *can* solve this challenge quite easily with a humble yet powerful HP-25, in a few minutes. Let's begin:

**Given the equation y' = x^2 + y^2, with y(0) = 0, find an accurate value for y(2).**

**Step 1:**

We need to use some algorithm and corresponding

program to find numeric solutions for 1st order differential equations. Any will do, the HP-25C Application Programs does include one (Euler's method), the HP-41C Advantage ROM includes a better one (4th order Runge-Kutta method), and there are may others available.

As I'm going to use my everyday-use HP-15C, this is the program I wrote from scratch, ad-hoc for this task. It's unoptimized, hastily written, but it does the job. It implements Runge-Kutta 4th order method, in 45 steps.

The 4th-order RK algorithm goes like this:

Given y' = f(x,y), with y0 = y(x0), we use a suitably

small 'step' size, h, to compute additional values, (x1, y1), (x2, y2), ...,

like this:
x1 = x0 + h

y1 = y0 + (k1 + 2k2 + 2k3 + k4)/6

where k1 = h.f(x0, y0)

k2 = h.f(x0 + h/2, y0 + k1/2)

k3 = h.f(x0 + h/2, y0 + k2/2)

k4 = h.f(x0 + h, y0 + k3)

once (x1, y1) have been computed, you then go on

likewise to compute (x2, y2), etc.

and the resulting HP-15C implementation is this:

01 LBL A 13 STO 3 25 RCL 1 37 GSB 2

02 STO 0 14 GSB 0 26 PSE 38 STO+ 3

03 Rdn 15 GSB 0 27 RCL 2 39 RTN

04 RTN 16 RCL+ 2 28 PSE 40 LBL 2

05 LBL B 17 RCL 0 29 GTO 1 41 RCL+ 1

06 STO 1 18 GSB 2 30 LBL 0 42 GSB C

07 X<>Y 19 RCL 3 31 2 43 RCLx 0

08 STO 2 20 6 32 / 44 STO+ 3

09 LBL 1 21 / 33 RCL+ 2 45 RTN

10 X<>Y 22 STO+ 2 34 RCL 0 46 LBL C

11 GSB C 23 RCL 0 35 2 47 RTN

12 RCLx 0 24 STO+ 1 36 /

where you define your function f(x,y) under LBL C,

assuming y is in the Y register and x is in the X register.

Set User mode and specify the step value h, like this:

h, A

Then enter the initial condition (x0,y0) and proceed to compute (x1,y1), (x2, y2), etc, like this:

y0, ENTER, x0, B -> (x1) -> (y1) -> (x2) -> y(2)

press R/S to stop the program when done. By the way, this program

*doesn't use any unique HP-15C features*, apart from recall arithmetic (RCLx 0) which you can easily simulate (RCL 0, x), so

*it will run* essentially unchanged on an HP-11C, HP-34C, and most any models having subroutine capability and labels. A version for the HP-25 and similar GSB-less machines is available as well.

**Step 2:**

Now, we need to define our equation, y' = x^2 + y^2. To this purpose, insert these steps after 46 LBL C

47 X^2

48 X<>Y

49 X^2

50 +

**Step 3:**
Let's try to find y(2). In User mode, FIX 7, using h = 0.1 ( .1, A), and x0 = 0, y0 = 0 (0, ENTER, 0, B), we let

the program run until it arrives at x = 2, duly noting the final result and several intermediate ones for x = 1, 1.7, 1.8, and 1.9. To check the accuracy of the result, we repeat the procedure a couple of times, halving the step size to h = 0.05, then to h = 0.025. This is what we get:

x y (h=0.1) y (h=0.05) y (h=0.025)

--------------------------------------------------

0 0 0 0

1 0.3502337 0.3502320 0.3502319

1.7 2.9727649 2.9727984 2.9727975

1.8 4.6866510 4.6880373 4.6881254

1.9 9.5161319 9.5622688 9.5666582

2 71.5789968 124.2866329 194.9352142

**Step 4:**
Now, instead of blindly accepting the last value for y(2), we must evaluate the results. From the table above, it's obvious that while y(1) = 0.3502319 is accurate to 7 decimals, y(1.8) is accurate to no more than 4 decimals, y(1.9) only to 2 decimals at best, and y(2) isn't accurate to any given digits. Further, the y-values seem to be more than doubling for a small linear increase (h) of the x-values, and this exponential growth should make us suspect a pole (singularity) at or near x = 2.

**Step 5:**

If a singularity is present nearby, reducing the step size h further to h = 0.001, 0.00001, or smaller, will be no use. The running times will be astronomic and the accumulation of rounding errors to reach x = 2 after 200,000 steps or more will make the final result inaccurate anyway. Do we reach for Mathematica ? No way, it's just that a different approach is called for.

As it's well known, HP vintage calcs were *"designed by engineers for engineers"*. Though many people now tend to take the easy ride and use pure brute force (i.e., Mathematica's N[ ..., 100]) instead of neuron-power, engineers of old were expected to have some familiarity with Calculus, at least elementary Calculus, plus a fine intuition for difficulties. On suspecting a pole, a singularity, any engineer worth his/her salt would think along these lines:

"Aw, gee, we've got some nasty singularity nearby, which is spoiling our algorithm's accuracy. Let's remove it out of existence ! How ? Easy ... let's make some change of variable, such as to turn a singularity-ridden function, like this unknown y(x), into some smooth, singularity-free one ... how about the arctan(x) function ? It will return a finite value, Pi/2 at most, for any argument x, even if it goes to infinity. Let's try !"

**Step 6:**

Lo and behold, let's make a change of variable:

[1] y = tan(z)

so then we have:

[2] z = arctan(y)

its derivative will be:

[3] z' = y'/(1 + y^2)

isolating y' gives:

[4] y' = z'.(1 + y^2) = z'.(1 + tan^2(z))

but the original equation says:

[5] y' = x^2 + y^2 = x^2 + tan^2(z)

so equating both expressions [4] and [5], we get:

[6] z'.(1 + tan^2(z)) = x^2 + tan^2(z)

which, after isolating z' gives:

[7] z' = (x^2 + tan^2(z)) / (1 + tan^2(z))

which is the new differential equation (in [x, z, z'], instead of the original [x, y, y']) we must now solve,

subject to the initial condition:

z(0) = arctan(y(0)) = arctan(0) = 0

and once we find z(2) using our HP-15C program, we'll simply have:

y(2) = tan(z(2))

**Step 7:**

The new differential equation is programmed under LBL C (after removing the original equation) like this:

47 X^2 52 1

48 X<>Y 53 LASTX

49 TAN 54 +

50 X^2 55 /

51 + 56 RTN

and now, in RAD mode, using h = 0.1, then h = 0.05, we get

x y (h=0.1) y (h=0.05)

------------------------------------

0 0 0

1 0.3368811 0.3368813

1.7 1.2463033 1.2463031

1.8 1.3606417 1.3606412

1.9 1.4666491 1.4666485

2.0 1.5676491 1.5676489

which clearly shows that the singularity has been effectively removed, so we've got:

z(2) = 1.5676489

and thus:

**y(2)** = tan(z(2)) **= 317.7225457**

which has

*nearly 7 digits correct*, greater accuracy being possible by simply using a smaller step size, h. Further, the exact x value corresponding to the pole near x=2 can be found accurately with our little HP calculator using a single iteration of Newton's method, like this:

Xpole = 2 - (z(2)-Pi/2)/z'(2)

where we already have z(2) and can compute z'(2) by simply using the differential equation [7]. So, calling z(2) = 1.5676489 = k, we then compute:

**Xpole** = 2 - ( k - Pi/2) / z'(2)

= 2 - ( k - Pi/2) / ((4+tan^2(k)) / (1 + tan^2(k)) )

**= 2.0031473**

which

*is accurate to all 8 digits shown*.

Thusly, we've succedded in finding both the value of y(2) and the singularity, to 7-8 digits accuracy, using just our small HP calculator and a little bit of the elementary knowledge and basic intuition any engineer should have in spades.

If any of you did find some of this useful or even learnt something, I think I will have proved my point and achieved my goal.

Best regards from V.