Well I don’t think there is an easy solution for exactly what you propose. plot essentially just takes 2 arrays of x and y coordinates and draws lines between them. How should it know which line to omit?

Some other ideas:

Instead of using two plot calls you can mark the jump with the value NaN. IIRC NaNs break the line. So your picewise function could return NaN exactly at the jump (or a small neighbourhood) and if your are plotting enoigh points such that one of them is inside that window, then it should work out.

Instead of plot/plot! you could use scatter/scatter! which just draws dots at the given location. If you play a bit with markersize and number of points you probably can get something reasonable as well.

Maybe there is a package for piecewise functions that can help you?

Yes you’ll need to play with some setting to make scatter look better. Did you try e.g. with color=:blue, markersize=0.25 or something like that? I am on mobile and can’t play with these settings myself.

Edit: You probably want to disable the border (black line) around the markers. You can try markerstrokewidth=0 for that.

It would probably make sense to use the information about the intervals that you have, by incorporating them into a type for which you can make a recipe which then makes gaps at the correct places. If you just have a blackbox function outputting numbers, you cannot really know whether two adjacent points are separated by a singularity or not, could always be continuous on a very small scale I guess?

As @jules mention, since this is just some blackbox you would need to provide it yourself with your own recipe, or break it into NaNs. scatter is probably okay too. Here’s one other approach you could take, identifying singularities by large slopes:

using Plots
function fix_singularities(y, x; slope_tol=100)
_x = similar(x, 0)
_y = similar(y, 0)
for i in firstindex(y):(lastindex(y)-1)
yᵢ, yᵢ₊₁ = y[i], y[i+1]
xᵢ, xᵢ₊₁ = x[i], x[i+1]
slope = (yᵢ₊₁ - yᵢ) / (xᵢ₊₁ - xᵢ)
if abs(slope) < slope_tol
push!(_x, x[i])
push!(_y, y[i])
else
push!(_x, x[i], (x[i] + x[i+1]) / 2)
push!(_y, y[i], NaN)
end
end
return _y, _x
end
g(x) = (x < 0.5) ? (x) : (x - 0.5)
x = 0.0:0.0001:1.0
y = g.(x)
_y, _x = fix_singularities(y, x)
plot(_x, _y)

Maybe this is robust enough for your use case, if you are willing to recompute the x and y like this (you can of course define your own method for providing a function and a grid still).

Random thought: Maybe Autodiff could identify true discontinuities? Just computing a local gradient does of course not contain this information. I am not too familiar with how all of these different packages work, but maybe one of the approaches could detect when there is a switch to a different code branch or just know due to code analysis when a jump occurs.

Yes but this quite different from plotting values from a blackbox. In principle you could define a type in Julia that representa a piecewise function and define a plotting recipe that handles the discontinuity.

The said functions in Mathematica is not the same thing as in Plots.jl. Mathematica knows it is defined piecewise (because of its symbolic nature) and does not assume continuity whereas, in Plots.jl, it is always numeric and it almost never knows it’s symbolic structure unless a recipe with a custom type defines piecewise function structure. Hence, like Matlab and several other plotting libraries, discontinuities in the lines of plot are inferred by NaN values.